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Abstract 

A major challenge to the successful full-scale development of modem aerospace 
systems is to address competing objectives such as improved performance, reduced costs, 
and enhanced safety. Accurate, high-fidelity models are typically time consuming and 
computationally expensive. Furthermore, informed decisions should be made with an 
understanding of the impact (global sensitivity) of the design variables on the different 
objectives. In this context, the so-called surrogate-based approach for analysis and 
optimization can play a very valuable role. The surrogates are constructed using data drawn 
from high-fidelity models, and provide fast approximations of the objectives and 
constraints at new design points, thereby making sensitivity and optimization studies 
feasible. This paper provides a comprehensive discussion of the fundamental issues that 
arise in surrogate-based analysis and optimization (SBAO), highlighting concepts, methods, 
techniques, as well as practical implications. The issues addressed include the selection of 
the loss function and regularization criteria for constructing the surrogates, design of 
experiments, surrogate selection and construction, sensitivity analysis, convergence, and 
optimization. The multi-objective optimal design of a liquid rocket injector is presented to 
highlight the state of the art and to help guide future efforts. 
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Corrected approximation 
Center of the radial basis function 
Total variance 

Partial variance of i th variable x,- 

.1 A 

m derivative of the function / 

Expected value of the quantity 

Expected value of the quantity considering all datasets 

Expected value of the square of the bias error at design point x 

Expected value of the variance at design point x 

True function to be modeled 

True response at design point x 

Vector of responses at sampled design points 

(Model) Response at i‘ h design point 

Model estimation 

Vector of radially symmetric functions 
i* radially symmetric function 

Matrix of radial basis functions for sampled design points 

Identity matrix 

Hadamard matrix 

Symmetric kernel function matrix 

Loss function 

Number of design variables 

Number of iterations in iterative fractional factorial design 
Number of terms in polynomial regression 
Number of radial basis functions 
Number of sampled design points 

Projection matrix 

Number of levels for each variable 

Correlation vector between the sampled design points and prediction design 
point 

Correlation matrix among the sampled design points 

Adjusted coefficient of multiple determination in polynomial regression 
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Sign of A* variable in m* iteration 

Strength of orthogonal array 
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Vector of coefficients of linear combinations of radially symmetric functions 
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i th design point 
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j th basis function for i th sample design point 
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i th component of the weight vector 

Vector of coefficients in polynomial regression 

Coefficients of the basis functions in polynomial regression 

Estimate of coefficient vector p in polynomial regression 

Index of the orthogonal array 

Radius of the radial basis function 

Model appraisal 

Probability density function 

Regularization parameter 

Mean of the responses at sampled design points 

Variance of the response at sample design points 

Adjusted root mean square error in polynomial regression 
Degree of correlation among data points along k? h direction 


1. Introduction 

Major risks in the successful full-scale development of aerospace systems are related to 
effectively addressing the competing needs of improving performance, reducing costs, and 
enhancing safety. Typically, the analysis of the components of such systems, such as air- 
breathing or rocket propulsion devices is expensive, hindering the search for optimal designs. 
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Surrogate-based analysis and optimization (SBAO) has been shown to be an effective 
approach for the design of computationally expensive models such as those found in 
aerospace systems, involving aerodynamics, structures, and propulsion, among other 
disciplines. Successful applications include the multidisciplinary optimal design of aerospike 
rocket nozzles, injectors and turbines for liquid rocket propulsion, and supersonic business 
aircrafts. 

For example, SBAO has been applied to high speed civil transport (Knill et al. [1]), airfoil 
shape optimization (Rai and Madavan [2]-[3], Madavan et al. [4]), diffuser shape optimization 
(Madsen et al. [5]), supersonic turbine (Papila et al. [6], Shyy et al. [7]-[8]), and injectors 
(Shyy et al. [7]-[8], Vaidyanathan et al. [9]-[l 0], Goel et al. [1 1]). 

This review provides a comprehensive discussion of the fundamental issues that arise in 
surrogate-based analysis and optimization (SBAO), emphasizing the concepts behind the 
methods and techniques, and practical implications of an integrated approach to SBAO for 
aerospace components and systems. 

The first part of the review is structured following the key steps in surrogate modeling 
(Figure 1): 

1. Design of Experiments (DOE). The design of experiment is the sampling plan in design 
variable space. The key question in this step is how we assess the goodness of such 
designs, considering the number of samples is severely limited by the computational 
expense of each sample. A discussion of the most prominent approaches related to DOE in 
SBAO is presented in Section III. 
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2. Numerical Simulations at Selected Locations. Here, the computationally expensive 
model is executed for all the values of the input variables in the DOE specified in the 
previous step. 

3. Construction of Surrogate Model. Two questions are of interest in this step: a) what 
surrogate model(s) should we use (model selection) and, b) how do we find the 
corresponding parameters (model identification)? A formal description of the problem of 
interest is discussed in Section II. A framework for the discussion and mathematical 
formulation of alternative surrogate-based modeling approaches is the subject of Section 
IV. 

4. Model Validation. This step has the purpose of establishing the predictive capabilities of 
the surrogate model away from the available data (generalization error). Section V 
discusses schemes for estimating the generalization error for model validation purposes. 

Then, surrogate-based analysis and optimization are discussed in Section VI and VII, 

respectively. A case study associated with the multi-objective optimal design of a liquid- 
rocket injector is used to illustrate the different issues that arise when conducting SB AO and 
particular methods and techniques selected. 

2. Overview of Surrogate Modeling 

With reference to Figure 2, surrogate modeling can be seen as a non-linear inverse problem 
for which one aims to determine a continuous function (/ ) of a set of design variables from a 
limited amount of available data (f). The available data f while deterministic in nature can 
represent exact evaluations of the function / or noisy observations and in general cannot carry 
sufficient information to uniquely identify / (multiple surrogates may be consistent with the 
available data as illustrated in Figure 3). Thus, surrogate modeling deals with the twin 
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A 

problems of: a) constructing a model / from the available data f (model estimation), and b) 
assessing the errors e attached to it (model appraisal). A general description of the anatomy 
of inverse problems can be found in Snieder [12]. 

Hence, using the surrogate modeling approach the prediction of the simulation-based 

A 

model output is formulated as f p (x) = f(x) + e(x ) . The prediction expected value and its 
variance V(f ) are illustrated in Figure 4, with 6 being a probability density function. Note 

that in Figure 3 it is assumed that the expected value of e(x) is zero. 

Different model estimation and model appraisal components of the prediction have been 
shown to be effective in the context of SBAO (see for example, Simpson et al. [13], 
McDonald et al. [14], Chung and Alonso [15], Jin et al. [16]), namely polynomial regression 
(PRG), Gaussian radial basis functions (GRF), and (ordinary) Kriging (KRG) as described by 
Sacks et al. [17]. Model estimation and appraisal components of these methods are presented 
in Section IV. 

A 

A good paradigm to illustrate how particular solutions ( / ) to the model estimation 
problem can be obtained is provided by regularization theory (see for example, Tikhonov and 
Arsenin [18], and Morozov [19]), which imposes additional constraints on the estimation. 

More precisely, / can be selected as the solution to the following Tikhonov regularization 
problem: 




D^l^dx 


( 1 ) 


where H is the family of surrogate models under consideration, L(x) is a loss or cost 


function used to quantify the so called empirical error (e.g., L(x) = x 2 ), X is a regularization 
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parameter, and D m ^ represents the value of the m-derivative of the proposed model at 

location x. Note that D m (/ represents a penalty term; for example, if m is selected equal to 2, 
it penalizes high local curvature. Hence, the first term enforces closeness to the data (goodness 
of fit), while the second term addresses the smoothness of the solution with X (a real positive 
number) establishing the tradeoff between the two. Increasing values of X provide smoother 
solutions. The purpose of the regularization parameter X is hence to help implement Occam's 
razor principle [20] which favors parsimony or simplicity in model construction. A good 
discussion on statistical regularization of inverse problems can be found in Tenorio [21] . 

The quadratic loss function (i.e., L 2 norm) is the most commonly used in part because it 
typically allows easy estimation of the parameters associated with the surrogate model; 
however, it is very sensitive to outliers. The linear (also called Laplace) loss function takes the 
absolute value of its argument (i.e., Lj norm); on the other hand, the Huber loss function is 
defined as quadratic for small values of its argument and linear otherwise. The so called s 
loss function has received considerable attention in the context of the support vector 
regression surrogate (Vapnik [22], Girosi [23]) and assigns an error equal to zero 
(interpolation) if the true and estimated values are within an e distance. Figure 5 illustrates 
the cited loss functions. 

Methods for identifying the regularization parameter X are typically based on 
generalization error estimates (e.g., cross validation); as a result, they will be discussed in the 
context of Validation (Section V). 

3. Design of Experiments (DOE) 

As stated in the Introduction, the design of experiment is the sampling plan in design 
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variable space and the key question in this step is how we assess the goodness of such designs. 
In this context, of particular interest are sampling plans that provide a unique value (in 
contrast to random values) for the input variables at each point in the input space, and are 
model-independent; that is, they can be efficiently used for fitting a variety of models. 

With reference to a) in most applications, where the assumed model is in doubt (see 
Sections II and IV), and the data is collected from deterministic computer simulations, the 
primary interest is minimizing the bias error because the numerical noise is small, and a DOE 
is selected accordingly. Brief descriptions for bias and variance components of the empirical 
error, and the corresponding expressions for the particular case of average error formulation 
and quadratic loss function are provided next. In the following discussion it is assumed that 
the function values (data set) have some noise or random component in them (e.g., due to 
numerical noise). 

A 

Bias : quantifies the extent to which the surrogate model outputs (i.e.,/(x)) differ from the 
true values (i.e., /(x) ) calculated as an average over all possible data sets D. 

A 

Variance-, measures the extent to which the surrogate model /(x) is sensitive to particular 
data set D. Each data set D corresponds to a random sample of the function of interest. 

For an average error formulation with a quadratic loss function the expressions for bias and 
variance are shown in Equations (2) and (3), respectively. In both expressions E ADS denotes 
the expected value considering all possible data sets. 

£^,(i) = {£^[/(i)]-/00} ! (2) 

— &ADS [/« -£«*[/(*) l] ! (2) 

There is a natural trade-off between bias and variance. A surrogate model that closely fits a 
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particular data set (lower bias) will tend to provide a larger variance. We can decrease the 
variance by smoothing the surrogate model but if the idea is taken too far then the bias error 
becomes significantly higher. In principle, we can reduce both the bias (can choose more 
complex models) and the variance (each model more heavily constrained by the data) by 
increasing the number of points provided the latter increases more rapidly than the model 
complexity. 

In practice, the number of points in the data set is severely limited (e.g., due to 
computational expense) and often during the construction of the surrogate model a balance 
between bias and variance errors is sought. This balance can be achieved, for example, by 
reducing the bias error while imposing penalties on the model complexity (e.g., Tikhonov 
regularization). 

When polynomial regression is used, and the true model is assumed to be a polynomial of a 
given order, there is good theory for obtaining minimum bias designs (e.g., Myers and 
Montgomery [24], Section 9.2) as well as some implementations in low dimensional spaces 
(e.g., Qu et al. [25]). For the more general case, the bias error can be reduced through a DOE 
that distributes the sample points uniformly in the design space (Box and Draper [26] and 
Sacks and Ylvisaker [27] as referenced in Tang [28]) but computational expense does not 
allow to conduct dense full factorial designs so, instead, some forms of fractional factorial 
designs (adequately chosen fraction of the full factorial design) are used. The uniformity 
property in designs is sought by (issue b), for example, maximizing the minimum distances 
among design points (Johnson et al. [29]), or by minimizing correlation measures among the 
sample data (Iman and Conover [30], Owen [31]). Practical implementation of these strategies 
includes orthogonal arrays (OA, e.g., Hedayat et al. [32]) and Latin Hypercube sampling 
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(LHS, e.g., McKay et al. [33]); the former produces uniform designs but can generate 
particular forms of point replications while the latter does not produce replicates but can lack 
uniformity. As a result, OA-based LHS (Tang [28], Ye [34]) and other optimal LHS schemes 
(Palmer and Tsui [35], Leary et al. [36]) have been proposed. These practical implementations 
are discussed next. 

3.1. Latin Hyper cube Sampling (LHS) 

Stratified sampling ensures that all portions of a given partition are sampled. LHS (McKay 
et al. [33]) is a stratified sampling approach with the restriction that each of the input variables 
(x*) has all portions of its distribution represented by input values. Following McKay et al.’s 
notation, a sample of size N s can be constructed by dividing the range of each input variable 
into N s strata of equal marginal probability l/N s and sampling once from each stratum. Let 
us denote this sample by x[ J) , j=l, ... N s ; k=l,.., K. The sample is made of components of 

each of the x*‘s matched at random. Figure 6 illustrates a LHS design. 

While LHS represents an improvement over unrestricted stratified sampling (Stein [37]) it 
can provide sampling plans with very different performance in terms of uniformity (issue b) 
measured by, for example, maximum minimum-distance among design points, or by 
correlation among the sample data. Figure 7 illustrates this shortcoming; the LHS plan in 
Figure 7(c) is significantly better than that in Figure 7(a). 

3.2. Orthogonal Arrays (OA) 

These arrays were introduced by C.R. Rao in the late 40’s and can be defined as follows. 
An OA of strength t is a matrix of N s rows and N dv columns with elements taken from a set 

of q symbols, such that in any N s xt submatrix each of the q‘ possible rows occurs the same 
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number y (index) of times. The number of rows (N s ) and columns ( N dv ) in the OA definition 

represent the number of samples and variables or factors under consideration respectively. The 
q symbols are related to the levels defined for the variables of interest, and the strength t is an 
indication of how many effects can be accounted for (to be discussed later in this section) 
with values typically between 2 and 4 for real-life applications. Such an array is denoted by 
OA(N s ,N dv ,q,t ). Note that by definition a LHS is an OA of strength one, OA (N s ,N dv ,N s ,l). 
There are two limitations on the use of OA for DOE: 

• Lack of flexibility. Given desired values for the number of rows, columns, levels, and 
strength, the OA may not exist. For a list of available orthogonal arrays, theory and 
applications, see, for example, Owen [38], Hedayat et al. [32] and references therein. 

• Point replicates: OA designs projected onto the subspace spanned by the effective factors 
(most influential design variables) can result in replication of points. This is undesirable for 
deterministic computer experiments where the bias of the proposed model is the main concern. 

3.3. Optimal LHS , OA-based LHS, Optimal OA-based LHS 

Different approaches have been proposed to overcome the potential lack of uniformity of 
LHS. Among those, most of them adjust the original LHS by optimizing a spreading measure 
(e.g., minimum distance or correlation) of the sample points. The resulting designs have been 
shown to be relatively insensitive to the optimal design criteria (Ye et al. [39]). Examples of 
this strategy can be found in the works of Johnson et al. [29], hnan and Conover [30] and 
Owen [31]. Tang [28] presents the construction of strength t OA-based LHS which stratify 
each ^-dimensional margin and shows that they offer a substantial improvement over standard 
LHS. Other strategies optimize a spreading measure of the sample points, but restrict the 
search of LHS designs which are orthogonal arrays, resulting in so called optimal OA-based 
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LHS (Leary et al. [36]). Adaptive DOE in which model appraisal information is used to place 
additional samples have also been proposed (Jones et al. [40], Sasena et al. [41], Williams et 
al. [42]). 

If feasible, two sets of DOE are generated, one (so called training data set) for the 
construction of the surrogate (Section rV) and one for assessing its quality (Validation as 
discussed in Section V). 

Given the choice of surrogate, the DOE can be optimized to suit a particular surrogate. This 
has been done extensively for minimizing variance in polynomial regression (e.g., Myers and 
Montgomery [24], Chapter 8) and to some extent for minimizing bias (e.g., Qu et al. [25]). 
However, for non-polynomial models, the cost of the optimization of a surrogate-specific 
DOE is usually prohibitive, and so is rarely attempted. 

4. Construction of the Surrogate Model 

There are both parametric (e.g., polynomial regression, Kriging) and non-parametric (e.g., 
projection-pursuit regression, radial basis functions) alternatives for constructing surrogate 
models. The parametric approaches presume the global functional form of the relationship 
between the response variable and the design variables is known, while the non-parametric 
ones use different types of simple, local models in different regions of the data to build up an 
overall model. This section discusses the estimation and appraisal components of the 
prediction of a sample of both parametric and non-parametric approaches, and the rationale of 
addressing the issue of model uncertainty (not knowing which surrogate may perform the 
best) using multiple surrogates (Zerpa et al. [43]). 

Specifically, the model estimation and appraisal components corresponding to Polynomial 
Regression (PRG), Kriging (KRG), and Radial Basis Functions (RBF) surrogate models are 
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discussed next, followed by a discussion of a more general non-parametric approach called 
Kernel-based regression. Throughout this section a square loss function is assumed unless 
otherwise specified, and given the stochastic nature of the surrogates, the available data is 
considered a sample of a population. 


4.1. Polynomial Regression model (PRG) 

The regression analysis is a methodology that studies the quantitative association between a 
function of interest f and N PRG basis functions zj, where there are N s sample values of the 

function of interest fi, for a set of basis functions z'° (Draper and Smith [44]). For each 
observation i a linear equation is formulated: 

M z) = "%/)/;>+ e, E(s,) = 0 V{s,)- a 2 (4) 

M 


where the errors are considered independents with expected value equal to zero and 
variance a 2 . For example, a second-order polynomial (N dv = 2; N PRG =6) can be expressed as: 


N AU 


/( x ) = fio + X fi * +XX Pij x i x j 


( 5 ) 


1=1 


i=l jii 


The estimated parameters p (by least squares) are unbiased and have minimum variance. 
The set of equations specified in Equation (4) can be expressed in matrix form as: 
f = Xp + s E(s) = 0 V(e) = a 2 I (6) 

where X is a N, x N PRG matrix of basis functions with the design variable values for 
sampled points. The vector of the estimated parameters then can then be found as: 

p = (X T X)‘'x T f (7) 
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Considering a new set of basis function vector z for design point P, the variance of the 


estimation is: 

(/(*)) -^(^(X'x)"'*) 


( 8 ) 


4.2. Kriging Modeling (KRG) 

It is named after the pioneering work of D.G. Krige (a South African mining engineer), and 
was formally developed by Matheron [45]. More recently. Sacks et al. [46]-[47], and Jones et 
al. [40] made it well-known in the context of the modeling, and optimization of deterministic 
functions, respectively. The Kriging method in its basic formulation estimates the value of a 
function (response) at some unsampled location as the sum of two components: the linear 
model (e.g., polynomial trend) and a systematic departure representing low (large scale) and 
high frequency (small scale) variation components, respectively. 

The systematic departure component represents the fluctuations around the trend, with the 
basic assumption being that these are correlated and the correlation depends only on the 
distance between the locations under consideration. More precisely, it is represented by a zero 
mean, second-order, stationary process (mean and variance constant with a correlation 
depending on a distance) as described by a correlation model. 

Hence, these models (Ordinary Kriging) suggest estimating deterministic functions as: 

f p (\) = ^ + e(\) E(e) = 0 cov(e(x (i) ),e(x u) ))*0 Vi,y (9) 

where fj. is the mean of the response at sampled design points, and e is the error with zero 
expected value, and with a correlation structure that is a function of a generalized distance 
between the sample data points. A possible correlation structure (Sacks et al. [46]) is given by: 


C0v(£'(x (i) ),£(x 0) )) = cr 2 exp 




) 


( 10 ) 
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where denotes the number of dimensions in the set of design variables x; cr, identifies 
the standard deviation of the response at sampled design points, and, <j> k is a parameter which 
is a measure of the degree of correlation among the data along the k th direction. Specifically, 
given a set of N s input/output pairs (x, f), the parameters, fx, a, and <f> are estimated such that 

a likelihood function is maximized (Sacks et al. [46]). Given a probability distribution and the 
corresponding parameters, the likelihood function is a measure of the probability of the 
sample data being drawn from it. The model estimates at unsampled points is: 

E(f p (x)) = {i + r T R~'(f-lju) (11) 

where the bar above the letters denotes estimates, r identifies the correlation vector 
between the set of prediction points and the points used to construct the model, R is the 
correlation matrix among the N s sample points, and 1 denotes an N s -vector of ones. 


On the other hand, the estimation variance at unsampled design points is given by: 


/ (l-l r R“‘r) 

V(f p (x)) = cr 2 l~r r R- r + - 1 


Gaussian Processes (Williams and Rasmussen [48], Gibbs [49]), another well-known 
approach to surrogate modeling, can be shown to provide identical expressions for the 
prediction and prediction variance to those provided by Kriging, under the stronger 
assumption that the available data (model responses) is a sample of a multivariate normal 
distribution (Rao [50]). 

4.3. Radial Basis Functions (RBF) 

Radial basis functions have been developed for the interpolation of scattered multivariate 
data. The method uses linear combinations of radially symmetric functions, hj(x), based 
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on Euclidean distance or other such metric, to approximate response functions as, 

jV /up- 

/„(*) = JwAW +f i ( 13 ) 

1=1 

where w represents the coefficients of the linear combinations, h/x) the radial basis 
functions and 5 independent errors with variance c? . 

The flexibility of the model, that is its ability to fit many different functions, derives from 
the freedom to choose different values for the weights. Radial basis functions are a special 
class of functions with their main feature being that their response decreases (or increases) 
monotonically with distance from a central point. The center, the distance scale, and the 
precise shape of the radial function are parameters of the model. 

A typical radial function is the Gaussian which, in the case of a scalar input, is expressed 
as: 

( x ~ c ) 2 

The parameters are its center c and its radius & . Note that the response of the Gaussian 
RBF decreases monotonically with the distance from the center, giving a significant response 
only in the center neighborhood. 

Given a set of N s input/output pairs (sample data) a radial basis functions model can be 
expressed in matrix form as, 

f = Hw (15) 

where H is a matrix given by, 


K (x) = exp 
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( 16 ) 


/t,(x (1) ) A,(x (,) ) ••• ^(x (,) ) 

R= A(x (2) ) ^(x (2) ) - ^(* (2) ) 

A(x W) ) ^(x W) ) - V,( xW ). 


Similarly to the polynomial regression method, by solving equation (15) the optimal 
weights (in the least squares sense) can be found to be, 

w = A-‘H T f (17) 

where A' 1 is a matrix given by, 

A* 1 = (H T H) * (18) 


The error variance estimate can be shown to be given by, 

f r P 2 f 


<j 2 = 


«r(P) 


where P is a projection matrix, 

p = i-ha‘h t 

The RBF model estimate for a new set input values is given by, 

/(x) = h r w 

where, h is a column vector with the radial basis functions evaluations, 


h = 


A,(x) 

h 2 (x) 

( X )J 


(19) 


( 20 ) 


( 21 ) 


( 22 ) 


A 

On the other hand, the prediction variance is the variance of the estimated model / (x) 


plus the error variance and is given by: 
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( 23 ) 


r(f, W) - >> r *)+ v (*) ■ =(h r (H r H)-‘ h+l)^ 

5 R 

The basic idea of RBF can be generalized to consider alternative loss functions and basis 
functions, in a scheme known as kernel-based regression. With reference to equation (1), it 
can be shown that independent of the form of the loss function L, the solution of the 
variational problem has the form (the Representer Theorem; see Girosi [23], Poggio and 
Smale [51]): 

/(x) = f> ( X(x,x (,) ) (24) 

i=l 

where K is a (symmetric) kernel function that determines the smoothness properties of the 
estimation scheme. Table 1 shows the kernel functions of selected estimation schemes with 
the kernel parameters being estimated by model selection approaches (see next section for 
details). If the loss function L is quadratic, the unknown coefficients in Equation (24) can be 
found by solving the linear system: 

(W J AI + K)a=f (25) 

where I is the identity matrix, and K a square positive definite matrix with elements 
K ij = K(x (0 ,x o) ) . Note that the linear system is well posed since (N 5 /II + K) is strictly 

positive and well conditioned for large values of N S A . If loss function L is non-quadratic, the 
solution of the variational problem still has the form of Equation (24) but the coefficients a, 

are found by solving a quadratic programming problem in what is known as support vector 
regression (Vapnik [22]). Comparative studies have shown that depending on the problem 
under consideration, a particular modeling scheme (e.g., polynomial regression, Kriging, 
radial basis functions) may outperform the others and in general, it is not known a priori 
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which one should be selected. See for example, the works of Friedman and Stuetzle [52], 
Yakowitz & Szidarovsky [53], Laslett [54], Giunta et al. [55], Simpson et al. [13], Jinet et al. 
[16]. Considering plausible alternative surrogate models can reasonably fit the available data, 
and the cost of constructing surrogates is small compared to the cost of the simulations, using 
multiple surrogates may offer advantages compared to the use of a single surrogate. Recently, 
multiple surrogate-based analysis and optimization approaches have been suggested by Zerpa 
et al. [43], based on the model averaging ideas of Perrone and Cooper [56], and Bishop [57]. 
The multiple surrogate-based analysis approach is based on the use of weighted average 
models which can be shown to reduce the prediction variance with respect to that of the 
individual surrogates. The idea of multiple surrogate-based optimization is discussed in a later 
section of the paper. 

5. Model Selection and Validation 

Generalization error estimates assess the quality of the surrogates for prediction and can be 
used for model selection among alternative models and establish whether they are ready to use 
in analysis and optimization studies (validation). This section discusses the most prominent 
approaches in the context of SBAO. 

5.1. Split Sample (SS) 

In this scheme the sample data is divided into training and test sets. The former is used for 
constructing the surrogate while the latter, if properly selected, allows computing an unbiased 
estimate of the generalization error. Its main disadvantages are that the generalization error 
estimate can exhibit a high variance (it may depend heavily on which points end up in the 
training and test sets) and that it limits the amount of data available for constructing the 
surrogates. 
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5.2. Cross Validation (CV) 

It is an improvement on the split sample scheme that allows the use of most if not all of the 
available data for constructing the surrogates. In general, the data is divided into k subsets (k- 
fold cross-validation) of approximately equal size. A surrogate model is constructed k times, 
each time leaving out one of the subsets from training, and using the omitted subset to 
compute the error measure of interest. The generalization error estimate is computed using 
the k error measures obtained (e.g., average). If k equals the sample size, this approach is 
called leave-one-out cross-validation (known also as PRESS in the polynomial regression 
terminology). Equation (26) represents a leave-one-out calculation when the generalization 
error is described by the mean square error (GMSE). 

GMSE = ]-£ l (f l - fry (26) 

K i=i 

where ff~ i} represents the prediction at x (,) using the surrogate constructed using all sample 

points except (x (,) , f t ). Analytical expressions are available for that case for the GMSE 

without actually performing the repeated construction of the surrogates for both polynomial 
regression (Myers and Montgomery [24], Section 2.7) and Kriging (Martin and Simpson [58]). 

The advantage of cross-validation is that it provides nearly unbiased estimation of the 
generalization error and the corresponding variance is reduced (when compared to split- 
sample) considering that every point gets to be in a test set once, and in a training set k-l times 
(regardless of how the data is divided); the variance of the estimation though may still be 
unacceptably high in particular for small data sets. The disadvantage is that it requires the 
construction of k surrogate models; this is alleviated by the increasing availability of surrogate 
modeling tools. A modified version of the CV approach called GCV-generalized cross 
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validation, which is invariant under orthogonal transformations of the data (unlike CV) is also 
available (Golub et al. [59]). 

If the Tikhonov regularization approach for regression is adopted, the regularization 
parameter X can be identified using one or more of the following alternative approaches: CV- 
cross validation (leave-one-out estimates), GCV (smoothed version of CV), or the L-curve 
(explained below). While CV and GCV can be computed very efficiently (Wahba [60], 
Hutchinson and de Hoog [61]), they may lead to very small values of X even for large 
samples (e.g., very flat GCV function). The L-curve (Hansen [62]) is claimed to be more 
robust and have the same good properties of GCV. The L-curve is a plot of the residual norm 


(first term) versus the norm of the solution for different values of the regularization 

J H 


parameter and displays the compromise in the minimization of these two quantities. The best 
regularization parameter is associated with a characteristic L-shaped “comer” of the graph. 


5.3. Bootstrapping 

This approach has been shown to work better than cross-validation in many cases (Efron 
[63]). In its simplest form, instead of splitting the data into subsets, subsamples of the data are 
considered. Each subsample is a random sample with replacement from the full sample, that is, 
it treats the data set as a population from which samples can be drawn. There are different 
variants of this approach (Hall [64], Efron and Tibshirani [65]) that can be used not only for 
model identification, but also for identifying confidence intervals for surrogate model outputs. 
This may come, though at the expense of considering several dozens or even hundreds of 
subsamples. 

For example, in the case of polynomial regression (given a model) regression parameters 
can be estimated for each of the subsamples and a probability distribution (and then 
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confidence intervals) for the parameters can be identified. Once the parameter distributions 
are estimated confidence intervals on model outputs of interest (e.g., mean) can also be 
obtained. 

Bootstraping has been shown to be effective in the context of neural network modeling; 
recently, its performance in the context of model identification in regression analysis is also 
being explored (Ohtani [66], Kleijnen and Deflandre [67]). 

6. Sensitivity Analysis 

Sensitivity, in this context, is a measure of the contribution of an independent variable to 
the total variance of the dependent data. Sensitivity analysis allows addressing settings such 
as: 

• Can we safely fix one or more of the input variables without significantly affecting the 
output variability (Screening)? 

• How can we rank a set of input variables according to their contribution to the output 
variability (Variables Prioritization)? 

• If we could eliminate the uncertainty of one or more of the input variables which ones 
should be chosen (Variable Selection for Maximum Uncertainty Reduction)? 

• If and which (group of) parameters interact with each other (Parameter Interactions)? 

• What are the main regions of interest in the parameter space if additional samples become 
available? 

• Does the model reproduce well known behavior of the process of interest (Model 
Validation)? 

There are alternative approaches for sensitivity analysis, differing, for example, in scope 
(local vs. global), nature (qualitative vs. quantitative), and in whether they assume a particular 
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model. Table 2 shows a sample of such methods. In this section, we discuss the Morris 
method (Morris [71]), Iterated Fractional Factorial Design (Andres and Hajas [72], Saltelli et 
al. [73]) frequently used for screening purposes and the Sobol’s method (Sobol [74]) a 
variance-based non-parametric approach for analysis applications. Note that a screening effort 
is expected to be economical (require a relatively small number of computationally expensive 
simulations) and potentially reduce the number of variables considered in the construction of 
the surrogate model. 

6.1. Morris Method (Morris [71]) 

This approach is model independent and is particularly useful when the number of 
available simulations for screening purposes is of the order of the number of design variables. 
Its main purpose is to determine, within reasonable uncertainty, whether the effect of 
particular design variables are negligible, linear and additive, or non-linear (interaction with 
other design variables are present). 

In its simplest form, the empirical distribution of the sensitivities associated with each of 
the design variables (x,) is estimated by computing a number ( 2N dv p ) of first order 

sensitivities of the model response with respect to each of the N dv design variables at a set of 

p random locations. Each x,- is then characterized by measures of central tendency (e.g., mean) 
and spread (e.g., standard deviation). A large (absolute) value of central tendency for x,- 
indicates a design variable with an important (global) influence on the model response. On the 
other hand, a large measure of spread is indicative of a design variable whose influence is 
highly dependent on the value of the other design variables (interactions/non-linear effects). 
While this formulation of the Morris method requires N s = 2 N dv p simulations, there are 
variations that result in more economical designs with some of the simulations used in 
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computing more than one sensitivity. A further reduction of the number of simulations can be 
obtained using cluster sampling (a concept to be discussed next). Details can be found in 
Morris [71]. 

6.2. Iterated Fractional Factorial Design (IFFD) Method (Andres and Hajas [72], 
Saltelli et al. [73]) 

The IFFD method was designed for screening purposes (identifying influential variables) 
for situations requiring fewer runs than there are design variables with the model response 
dominated by a few influential variables. Assuming a particular model, IFFD establishes as 
influential design variables with significant linear, quadratic, or non- linear/interaction effects. 

IFFD belongs to the so called cluster sampling designs with the design variables randomly 
aggregated into a small number of clusters. These clusters are then investigated using 
orthogonal fractional factorial designs in multiple iterations (composite design), with different 
groupings of design variables. A fractional factorial design is defined as a factorial experiment 
in which only an adequately chosen fraction of the combinations required for the complete 
factorial experiment is selected to be run while the orthogonal property makes reference to 
balanced designs in which the number of different combination of values for two or three 
design variables appear with equal frequency. Considering an influential group must contain 
an influential design variable, the influential variables are expected to lie in the intersection of 
influential groups in the iterations. 

The steps involved in the IFFD method can be summarized as follows: 

1. The design variables are sampled at three different levels, L (low), M (medium) and H 
(high) ensuring that the sampling is balanced. The number of clusters or groups of 
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variables k is set. The value of k is assumed a power of two without a loss of generality; 
commonly used values of k are 8 and 16. 

2. A basic two-level (L&H) resolution IV design matrix is constructed. The resolution of a 
design makes reference to its ability to discriminate between main factors (the individual 
effect of design variables), two-factor interactions, three factor interactions, etc. In a 
resolution IV design, main factors are confounded at worst with three-factor interactions; 
two factor interactions are confounded with certain other two-factor interactions though. A 
Hadamard matrix can describe these designs with the value of -1 or +1 at each entry 
representing the L and H values, respectively. Each column in the matrix represents the 
values of the variables assigned to it, while each row denotes a computational experiment 
to be conducted. This Hadamard matrix will be denoted by Jk[i, j] and has a size of 2k 
rows and k columns. 

3. Randomly assign (with equal probability) a design variable to a column in the two-level 
resolution IV design matrix; the corresponding variable value will be given by the design 
matrix multiplied by a sign (i.e., +1 or -1 with equal probability). The sign is denoted by 

s ” ; the sign of variable A in iteration m. 

4. Repeat step 3 N ilr times assuming the number of available simulations for screening 
purposes is approximately 2 kN ilr . During a fraction of the iterations (usually !4) each of 

the design variables values are set to zero, hence creating a three-level design, even though 
each of the iterations has either a single middle level M or two extreme levels (L, H) per 
individual variable. 
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5. The model responses corresponding to the rows in the design matrix in each of the 
iterations are calculated. The model response in iteration m for a particular row is denoted 


by fr- 

The global effect of a variable Xj that takes its value from column j in the design matrix in a 
given iteration can be calculated as: 




(27) 


The global effect of a variable denoted with letter A throughout the entire design is given 


by: 


fXMEtv.f') 


ME(A,t) = ,f) | jJ * 0) = 


m=l 


M 


m = 1 


(28) 


where C* is the column that has assigned variable A in iteration m. Quadratic effects are 


estimated by: 

QE(A,f) = avg( f 1^=0)- avg{ f \s*0) 

N* , , Ik N* 2k 


QE(A,() = 


m = 1 


i=l 


(29) 


itr . n >tr 


m= 1 


m=l 


Provided the model response is dominated by a few influential variables, and their 
contribution mostly limited to linear and quadratic effects, the IFFD method has been found to 
be effective in the screening of hundred and even thousands of variables with an order of 
magnitude less simulations. 


6.3. Sobol’s Method (Sobol [74]) 

To understand the concept, assume a surrogate model of a square integrable objective, /(x), 
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as a function of a vector of design variables, x, whose values have been scaled between zero 
and one (this assumes that the design domain is box-like). This surrogate model represented in 
Equation (30) can be decomposed as the sum of functions of increasing dimensionality as: 

/ M = fo + Z ft fa ) ' + E fv (*< » ■ x j ) ' + - ' + fn...K (*i , * 2 » •••> x K ) (30) 

i i<j 

If the following condition 

JA,A =0 < 31 > 

o 

is imposed for k = //, i s , where 1 < i x < ... < i s < N dv , the decomposition described in 

Equation (30) is unique and each term in the sum can be obtained by computing the following 
integrals: 


j/wn 

(32) 

|/wn*i =fo + fi( x i) 

k*i 

(33) 

from which /{x,) can be found, and 


\f ( x ) I! ^ = fo + fi ( x - ) ' + fj ( x j ) ' + fu ( x < ’ X J ) 

(34) 


k*i,j 


from which fi/xu xj) can be obtained. The higher dimensional summands are similarly 
foxmd except for the last one that is calculated using Equation (30). 

Furthermore, the summands are orthogonal (ensured by the condition expressed in 
Equation (31)), and, square integrable following the same condition for j{x). Therefore the 
partial variances, that is, the contribution of each of the summands to the total variance 
observed in the response, can be shown to be: 


D. . = f f 2 , dx, ...dx, 

it-', j J «)•••'» <i *. 


(35) 
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with the total variance being equal to: 


D= jf 2 (x)dx-f 2 

which can also be expressed as, 


AU N* 


*>-! I A l 4 

5=1 ij 


( 36 ) 


(37) 


Each partial variance gives a measure of the contribution of each independent variable or 
set of variables to the total variance, and provides an indication of their relative importance. 
Note that all the required integrations are conducted on the surrogate (fast) model and can in 
principle be calculated accurately provided an integration numerical procedure is available 
(e.g., Gaussian quadrature). 

The relative importance of a design variable is quantified by a set of indices, namely, main 
factors, and total sensitivity indices. The former refer to the fraction of the total variance 
contributed by a particular variable in isolation, while the latter represents the contribution 
(expressed as a fraction) of all the partial variances in which the variable of interest is 
involved. The influence of a design variable, say x h to an objective variability without 
accounting for any of its interactions with other variables is denoted as a main factor index 
and given as: 

S i =D i /D (38) 

To calculate the total sensitivity of any design variable, say x h the design variable vector x 
is divided into two complementary subsets, and Z where Z is a vector containing X\,X 2 , X 3 , 
..., x„ (n ^ i). The purpose of using these subsets is to isolate the influence of Xj on the J[x) 
variability from the influence of the remaining design variables included in Z. The total 
sensitivity index for Xi is then defined as 
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(39) 


g tota i jy ota ^ I d 

where 

D™ 1 =D i+ D iZ (40) 

£> is the partial variance associated with x iy and D iZ is the sum of all partial variances 

associated with any combination of the remaining variables representing the interactions 
between *,• and Z. Similarly the sum of all the partial variances associated with the variables 
denoted by Z can be defined as Dz, hence the total variance can be written as, 

D = D i+ D z +D. z (41) 

Formulations of the Sobol's method that account for non-rectangular domains and 
correlated inputs are available; see, for example, Jacques et al. [75], and Mack et al. [76] for a 
recent application. A detailed discussion of global sensitivity methods and applications can be 
found in Sobol [74], Homma and Saltelli [77], Saltelli et al.[78] and the references therein. 

7. Surrogate-based Optimization 

Surrogate model based optimization refers to the idea of speeding optimization processes 
by using surrogates for the objectives and constraints functions. The surrogates also allow for 
the optimization of problems with non-smooth or noisy responses, and can provide insight 
into the nature of the design space (Vaidyanathan et al. [10], Goel et al. [1 1], Mack et al. [76], 
Goel et al. [79]). SBAO has shown to be effective in both multidisciplinary and multi- 
objective optimization. This section discusses the basic unconstrained SBAO algorithm, a 
newly proposed multiple surrogate-based optimization approach, the use of surrogate 
management frameworks to obtain provably convergent methods, and ways of addressing 
general non-linear constraints. 


30 



7.1. Basic Unconstrained SB AO 


It can be summarized as follows: 

1 . Construct a surrogate model from a set of known data points 

2. Estimate the function minimizer using the surrogate function 

3. Evaluate the true function value at the estimated minimum ( checking phase ) 

4. Check for convergence; if achieved, stop. 

5. Update surrogate using new data points 

6. Iterate until convergence 

where possible convergence criteria include the achievement of a target for the objective 
function, or finding that the best point found so far remains unchanged (within a given 
tolerance) from one iteration to the next. There are trust-region approaches, to be discussed 
below where the region where the surrogate is constructed and the optimization is conducted, 
is shrunk or expanded depending on the results of the checking phase in order to guarantee 
convergence to a local optimum (e.g., Alexandrov [80]). In the so called one-shot solution 
approach though, a single iteration is used. 

In general, the algorithm above should be focused on identifying trends in the objective 
with the accuracy of the surrogates being important only when in the vicinity of a minimizer. 
Estimating the function minimizer is conducted using standard optimization methods and 
updating the surrogate using new data points can be accomplish by, for example, using model 
appraisal information and merit functions (e.g., Jones et al. [40], Torczon and Trosset [81], 
Sasena et al. [41], Queipo et al. [82]-[84]). Jones et al. [40], for example, use a generalized 
expected improvement function whereby points with low objective function value or high 
uncertainty are given precedence. A discussion of other sampling criteria in the context of 
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SBAO is presented by Sasena et al. [41 ] as proposed by Watson and Barnes [85]. 

7.2. Multiple Surrogates SBAO 

For most challenging problems, given the sparseness of the data used to construct the 
surrogates, alternative surrogates can provide reasonable approximations while giving 
different uncertainty estimates throughout the design space. On the other hand, the surrogate 
construction requires small computational resources compared to the cost of simulations and, 
as previously discussed, a variety of parametric and non-parametric alternative loss functions 
(e.g., quadratic, and e -sensitive) can be implemented. 

The use of multiple surrogates will bring in several advantages at the optimization and 
decision-making levels of the proposed approach, namely: 

• The predictions made by each surrogate are checked in the checking phase and when 
simulations are performed in additional cycles (iterations). Thus, we should be able to rank 
surrogates based on a small number of cycles and select those that appear to fit best the 
problem at hand. 

• The additional surrogate optima analyzed in the checking phase may increase the global 
search capability of the optimization in that they may correspond to different local optima. 

• A weighted averaged model may be constructed that is expected to provide a prediction 
with lower variance than any of the individual surrogates. 

• Large variability in the estimated values and variances among surrogates at a point in 
design space may indicate that the uncertainty at that point is higher than predicted. 

Some preliminary results obtained on the use of multiple surrogates for optimization has 

been reported by Zerpa et al. [43]. 

With reference to the basic SBAO algorithm, when multiple surrogates are available 
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changes need to be made in the proposed approach. At the optimization level, in general there 
will be multiple suggested optima, and the checking phase will include not one but as many 
simulations as surrogates. This new data set is available to all surrogates and their expressions 
for the predictions and variances can be updated. At the decision making level, the designer 
can make a qualitative assessment of the uncertainty in the results, and has more information 
to decide whether or not to undertake another cycle and on how to proceed (the number and 
location of additional sample points). Note that significant differences in the predictions and 
variances at particular locations can be used to identify regions where the uncertainty may be 
higher than expected. Hence, it may be required to sample those regions and subsequently 
update the predictions and variances. 

7.3. Convergence and SBAO 

There are two approaches that are provably globally convergent to local optima or 
solutions of the original (computationally expensive) optimization problem of interest that 
involve surrogates. These are called the Approximation Model Management Framework 
(AMMF) (Alexandrov [80], Lewis [86], Rodriguez et al. [87]) and Surrogate Management 
Framework (SMF) (Booker et al. [88]), and inherit/exploit the convergence properties of trust 
region and pattern search methods, respectively. They both use variable-fidelity models 
( surrogate and computationally expensive) and differ in whether they can use derivative 
information of the high fidelity function (if available). A description of each of these 
frameworks is next. 

7.4. Approximation Model Management Framework (AMMF) 

The approach is typically associated with gradient-based optimization algorithms and relies 
under the general assumption that the surrogate model is accurate enough for the purpose of 
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finding a good direction of improvement for the higher fidelity (computationally expensive) 
model. The AMMF replaces the local, Taylor series approximation, typical of conventional 
optimization, by an arbitrary model required to satisfy consistency conditions. The 
consistency conditions can be of different orders (zero, first, or second); for a first order model 
the consistency conditions are: 

a c (X) = /(x c ) 

Va c (xJ = V/(x e ) (42) 

where a and / denote the corrected approximation (its expression is shown below) and high 

fidelity models, respectively. The consistency conditions can be imposed using a correction 
technique: 

^(x) = A( x )/yi 0 (x) (43) 

where / w (x) and / /o (x) corresponds to computationally expensive and surrogate 

evaluations, respectively. Furthermore, a first order model p c of P about a current design 
variable vector x c can be written as: 

j3 c (x) = /3(x c ) + VJ3(x c ')(x-x c ) (44) 

The local model of P is then used to correct f, 0 (x) to obtain a better approximation a(x) of 

/«(*): 

/*(x) = y^(x)y; 0 (x) 

a(x)=A(x)y; 0 (x) (45) 

It can be shown that the approximation model a(x) satisfies the first order consistency 
conditions. The basic AMMF algorithm can then be summarized as follows: 

Until convergence do: 
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• Choose the corrected approximation a k so that it satisfies the consistency conditions 

• Find an approximate solution s k to the subproblem: 
min a k (x k + s) 

subject^ L k ( J 

• Compare the actual and predicted decrease in f u 

• Update x k and A k accordingly 

As in the classical trust region methods (Dennis and Schnabel [89]) the length of the steps 
is regulated based on how well the corrected model predicts the decrease in f hi . A significant 

mismatch in the actual and predicted decrease in f hi may lead to a reduction in the trust radius 

and/or to update the surrogate model approximation. The AMMF approach inherits 
convergence properties of the classical tmst regions algorithms for non-linear optimization. 
More precisely, first order AMMF methods can be shown to converge to local optima of the 
high fidelity problem under appropriate standard conditions of continuity and boundedness of 
the functions and derivatives, given that the consistency conditions are imposed at each major 
iteration. The AMMF can be adapted to handle nonlinear constraints and successful 
applications related to variable fidelity models have been reported (Alexandrov et al. [90], 
[91]). While the guaranteed convergence of the algorithm is valuable, it comes at the cost of 
requiring derivative calculation for the expensive simulation at one point in the design space 
(normally the optimum of the previous cycle). 

7.5. Surrogate Management Framework (SMF) 

It is a framework based on pattern search methods (see, for example, Torczon [92]) which 
incorporates surrogates to make the optimization cost effective. Being based on pattern search 
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methods, the SMF framework: i) can be used for scenarios where the objective functions are 
non-differentiable, or for which sensitivities are difficult or expensive to attain, ii) can be 
easily adapted to either a parallel or distributed computing environment, iii) is less likely to be 
trapped by non-global minimizers than are traditional nonlinear optimization algorithms, iv) 
can be extended to handle general nonlinear constraints, and v) convergence theory for the 
SMF can be found in the convergence of pattern search methods. The following SMF 
description follows that provided by Booker et al. [88] and Marsden et al. [93]. 

The SMF algorithm (illustrated in Figure 8) is mesh-based (all points evaluated are 
restricted to lie on a mesh) and consists of two steps, SEARCH and POLL. The exploratory 
SEARCH step uses the surrogate to aid in the selection of points likely to minimize the 
objective function. The SEARCH step provides means for local and/or global exploration of 
the parameter space, but it is not strictly required for convergence. Hence, this step can be 
adapted to the particular engineering problem under consideration. 

The convergence of the SMF is guaranteed by the POLL step, in which points neighboring 
the current best points on the mesh are evaluated in a positive spanning set of directions 
(positive basis) to check whether the current best point is a mesh local optimizer. A positive 
basis in is a set of vectors whose nonnegative linear combinations span R N * , but for 
which no proper subset has that property. The relevance of a positive basis in the context of 
the SMF is that it ensures that if the gradient of the high fidelity function / at the best current 
solution x is not zero, then at least one vector in the positive basis defines a descent direction 
from / at x. This can be guaranteed without any knowledge of the gradient. Any positive basis 
has at least Nd V +l (minimal) and at most 2N<j v (maximal) vectors; for example, in three 
dimensions, such a basis can be given by: (1,0,0), (0,1,0), (0,0,1), and (-1,-1,-!). For 
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unconstrained problems, a minimal positive basis is sufficient to guarantee convergence. The 
POLL set is made by the set of mesh points adjacent to the current best point in the above- 
referenced directions. 

In the SEARCH step, the high fidelity model is evaluated at one or more minimizes as 
predicted by the surrogate. If a lower objective function value is found among the recent 
evaluations, the search is considered successful, the surrogate is updated, and another search 
step is performed. If the SEARCH step fails to find an improved point, then it is considered 
unsuccessful and a POLL step is performed. In this step, the set of POLL points are evaluated. 
As soon as an improved point is found, a SEARCH step can be conducted on the current mesh. 
If no improved points are found, then the current best point is considered a so called mesh 
local optimizer. The local term acknowledges the fact that the search was limited to the POLL 
set. For greater accuracy the mesh can be refined, at which point the algorithm continues 
with a SEARCH. Convergence is reached when a local minimizer on the mesh is found, and 
the mesh has been refined to the desired accuracy. Each time new data points are found in a 
SEARCH or POLL step, the surrogate model is updated. 

Convergence of the SMF framework for constrained problems can be guaranteed (provided 
the function is continuously differentiable in the feasible region) if a generalized pattern 
search (GPS) strategy is adopted. The GPS strategy step halves the mesh if the POLL step is 
unsuccessful. Convergence properties of the SMF based on the generalized pattern search 
strategy can be found in Booker et al. [88]. Note, however, that the number of high-fidelity 
simulations performed (the POLL stage) is relatively small compared to those used for 
constructing surrogates, so that more iterations are required, and parallel processing 
opportunities are more limited. 
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7. 6. Constrained SB A O 


Once the constraints are modelled using surrogates, constrained SBAO is typically solved 
using well-known non-linear programming algorithms (Luenberger [94]) for constrained 
optimization based on the concept of penalty function -penalty method, classical and modified 
barrier methods, augmented Lagrangian methods- which approximate the constrained problem 
to an unconstrained one. The penalty function is a linear combination of the objective and 
some measure of the constraint violation. The approximation is accomplished in the case of 
penalty methods by adding to the objective function a term that prescribe a high cost for 
constraint violations and in the case of barrier methods by adding a term that favors points 
interior to the feasible region over those near the boundary. A related idea is that 
corresponding to an augmented Lagrangian function in which a penalty term is added to a 
Lagrangian function. While successful applications of these methods in the context of SBAO 
have been reported, in general, there are some difficulties associated with the use of penalty 
functions and the selection of the penalty parameter. Namely, there is usually a threshold 
value for the penalty parameter below which the unconstrained problem associated with the 
penalty function does not have a local minimum at the solution of the constrained problem, 
and the threshold value is not known a priori. If the penalty parameter is too low the solution 
can be unfeasible, and if it is too high it can damped out the effect of the objective which can 
result in very slow convergence. 

Methods to avoid the above-referenced difficulties have been proposed. In particular, 
Fletcher and Leyffer [95] recently presented a so called filter approach that does not require 
estimating a potentially troublesome penalty or barrier parameter. In this method, instead of 
combining the objectives and a measure of constraint satisfaction into a single minimization 
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problem (as with penalty functions), these are treated as two separate objectives (multi- 
objective optimization). The set of points that are non-dominated (none is better than any 
other in both objective function value and constraint violation) is called a filter , and an 
iteration of an algorithm, selects a filter point, generates a new one (e.g., using gradient 
information) and accept the point (added to the filter) if the new point is non-dominated with 
respect to those in the filter, or reject it otherwise. Note the contrast with the penalty function 
methods where a point in an iteration is accepted only if there is a decrease in the penalty 
function. While the filter approach was developed in the context of gradient-based 
optimization (e.g., SLP, SQP) the method has been extended to derivative-free methods 
(Audet and Dennis [96]) and holds promise to be useful in industrial applications (see, for 
example, Audet et al. [97]; Marsden et al. [98]). 

8. Case Study: Multi-objective Liquid-Rocket Injector Design 

In this section, we will use an injector design case study, motivated by rocket propulsion, to 
illustrate the issues and usefulness of SBAO for multi-objective optimization. The materials 
are extracted largely from the recent conference papers by Vaidyanathan et al. [9]-[10] and 
Goel et al. [11]. 

Two types of injectors are commonly used in space propulsion. These are (i) coaxial 
injectors and (ii) impinging type injectors. In coaxial injectors, the propellant streams flow in 
parallel. There are different types of coaxial injectors namely, shear coaxial, and swirl coaxial 
injectors. In the shear coaxial injector element mixing is accomplished through a shear-mixing 
process [99]. In impinging type injectors, the mixing occurs by direct impingement of the 
propellant streams at an acute angle. Calhoon et al. [99] conducted a large number of cold- 
flow and hot-fire tests over a range of propellant mixture ratios, propellant velocity ratios and 
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chamber pressure for shear coaxial, swirl coaxial, impinging, and premixed elements. The 
data were correlated directly with injector/chamber design parameters, which were recognized 
from both theoretical and empirical standpoints as the controlling variables. A schematic 
diagram of a single element injector proposed by Boeing is shown in Figure 9. 

8.1. Problem Description 

The injector design has two primary objectives: improvement of performance and life. As 
discussed by Vaidyanathan et al. [9] the performance of the injector is indicated by the axial 
length of the thrust chamber, while the life of the injector is associated with the thermal field 
inside the thrust chamber. A visual representation of the objectives is shown in Figure 10. 

In summary, the objectives are listed as the following: 

Combustion Length (A cc ): This is defined as the distance from the inlet where 99% of the 
combustion is complete. It is desirable to keep the combustion length as small as possible as 
this directly affects the size of the combustor and hence the weight of the spacecraft. 

Wall Temperature {TW 4 ): This is defined as the wall temperature at 3" from the injector 
face. Higher values of the wall temperature have negative impact on the life of the injector, so 
this objective has to be minimized. 

Face temperature ( TF max ): This is defined as the maximum temperature of the injector 
face. It is desirable to reduce temperature to increase the life of the injector. 

Tip temperature (TT max ): This is defined as the maximum temperature on the post tip of 
the injector. It is desirable to keep this temperature as low as possible. 

The independent variable ranges considered in this study are shown in Table 3. 

It can be seen that the dual goal of maximizing the performance and the life is now cast as 
a four-objective design problem. Improved performance of the injector causes development of 
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higher temperatures in the thrust chamber, so as discussed by Goel et al. [11], some of these 
objectives are conflicting in nature, hence there cannot be a single optimal solution for this 
problem. 

There are four primary design variables for the injector design problem shown in Figure 
10. These variables are given as follows: 

Flow Angle (a): This variable represents the hydrogen flow angle. The maximum angle 
varies between (f to ct. 

Hydrogen area {AHA): This variable represents the increment with respect to the baseline 
cross-section area of the tube carrying hydrogen. The increment varies from 0-25% of the 
baseline hydrogen area. 

Oxygen area (AOA): This variable represents the decrement with respect to the baseline 
cross-section area of the tube carrying oxygen. The area varies between 0 - (-40) % of the 
baseline area. 

Oxidizer post tip thickness ( OPTT ): Oxidizer post tip thickness varies between X" to 
2X", where X is a dimension, which cannot be disclosed because of confidentiality 
restrictions. 

All the variables are linearly normalized between 0 and 1. More information about the 
design variables can be found in Vaidyanathan et al. [9] 

The objective functions values corresponding to a particular design can be obtained 
through computationally expensive CFD simulations. Each CFD simulation can be obtained 
from a pressure-based Navier-Stokes solver, FDNS500-CVS [100]-[102]. In addition to the 
Favre-averaged Navier-Stokes equations, the two-equation turbulence model, and kinetic 
equations are solved. Steady state solutions are reached by an implicit Euler time marching 
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scheme. The chemical species transport equations represented the H 2 - O 2 chemistry with the 
aid of a 7-species and 9-reaction set [100]-[102]. The simulation domain and the boundary 
conditions used in all the CFD cases are shown in Figure 11. Because of the very large aspect 
ratio, both the injector and chamber have been shortened (at the cross hatched areas) for 
clarity. This is a computationally expensive simulation-based problem so for optimization 
purposes it is advisable to develop the surrogate models for the objective functions. It was 
shown by Vaidyanathan et al. [9], [103], Shyy et al. [8], and Papila et al. [6], [104] that 
accurate response surfaces for complex problems, like single element injectors, can be 
developed. 

Comparing two of the evaluated designs (using CFD), Vaidyanathan et al. [9] have shown 
the motivation for using surrogate models and an efficient optimization technique in the 
design process. The independent design variables, normalized between 0 and 1, are shown for 
the two cases in Table 4. In terms of the design space evaluated, these two designs (X and Y) 
are seen to be quite different. 

The chamber wall and injector face temperatures (TW 4 and TF max , respectively) for Case Y 
are low or moderately low, while for Case X, they are high. It was seen that a large 
recirculation zone located between the injector and the chamber wall strips hot gases from the 
flame and causes them to flow back along the chamber wall and injector face [9]. This 
phenomenon regulates the chamber wall temperature TW 4 and the injector face temperature 
TFmax. The other life-indicating variable, the maximum oxidizer post tip temperature, TT max 
has essentially the opposite trend as compared to the other two temperatures. The performance 
indicator, combustion length X cc is seen (Table 5) to be at a minimum level for Case X 
(shorter combustion lengths indicate better mixing elements) and at a moderate level for Case 
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Y. It is clear that the performance and environment indicators exhibit competing trends such 
that no design is the “best” for all the indicators of performance and environment. 

These comparisons confirm the earlier statement that changes in the injector design details 
have major effects on injector performance and injector-generated environments. These, and 
other relevant geometric details, can be effectively addressed by modem computation-based 
design tool. 

8.2. Design of Experiments (DOE) 

The selected DOE scheme was orthogonal arrays (Section III). Specifically, an OA (54, 4, 
3, 2) was used to generate fifty four (54) designs, considering four (4) factors, three (3) levels, 
with strength 2. For model selection and validation the split-sample approach (Section V) was 
selected so from the fifty four (54) designs, fourteen (14) were set aside for testing purposes. 
The fifty four (54) designs specified by the orthogonal array were computed on an 
axisymmetric geometry with 336x81 nodes. Only 33 out of the 40 training cases gave valid 
results. Results of the remaining seven cases contained unsteady features, which did not 
represent solutions of the steady-state model employed. Based on the quality of the 
approximation obtained using the selected surrogate model (polynomial regression) for TT max 
and X cc , it was noticed that the grid distribution in the combustion zone was insufficient. The 
grid was then refined to a 430x81 grid (Figure 12). The refined grid was found to be 
appropriate and was selected for the purposes of the optimization study. Only 2 out of the 
forty (40) designs to be used for fitting the surrogate model gave unacceptable results (i.e., 
unsteady features were present). Note that in addition to facilitating design optimization, 
surrogate models can also help check the adequacy of CFD simulations via identification of 
outliers. 
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Hence, the final data set included thirty eight (38) designs for fitting the surrogate model 
and fourteen (14) to test their predictive capabilities. The fitting and testing design points can 
be found in Vaidyanathan et al. [9]. 

8.3. Construction of the Surrogate Model 

The parametric polynomial regression approach (Section IV) was selected as surrogate 
model. Before constructing the surrogates, all the design variables were scaled between zero 
and one based on their upper and lower bounds. For example, when the AOA is 1 or 0, the O 2 
flow area is reduced by 0% or 40%, respectively, as compared to the baseline area. Similarly, 
the objective values were scaled between zero and one based on the upper and lower bounds 
observed in the sample data. Once the polynomial approximations were available, the scaled 
objective values were normalized using the surrogate-based minimum and maximum values 
for the objectives. 

The polynomial-based approximations were created using standard least-squares regression 
[24], and terms with insignificant influence (based on t-statistics) on the prediction of the 
response/objective were discarded, thereby improving the parsimony of the surrogate model. 
The statistical analysis software JMP [106] was used for the generation of the polynomials. 
The quality of alternative polynomial approximations was evaluated by comparing their 
adjusted root mean square (rms)-error <y a , and adjusted coefficient of multiple determinations, 
R a 2 [24] corresponding to the training set, and those with the best performance were selected 
(model selection). The model evaluation process was conducted using the rms-error a for the 
test set. Details of the model selection and validation (Section V) process can be found in 
Vaidyanathan et al. [9]. 

The polynomial approximations corresponding to each of the objectives are presented 
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below. Note that the best approximations for objectives TF max , TW 4 and X cc are quadratic 
polynomials, while for objective TT max a reduced cubic polynomial represented the best choice. 

TF max = 0.692 + 0.477(a) - 0.687(AHA) - 0.080(AOA) - 0.0650(OPTT) - 0.167(a ) 2 - 
0. 01 29 (AH A )(a) + 0.0796(AHA) 2 - 0.0634(AOA)(a) - 0.0257(AOA)(AHA) + 0.0877 (AOA) 2 - 
0.052 l(OPTT)(a) + 0.00 156(OPTT)( AHA) + 0.00 198(OPTT)( AOA) + 0.0184(OPTT) 2 (47) 

TW 4 = 0.758 + 0.358(a) - 0.807 (AHA) + 0.092 5 (AOA) - 0.0468(OPTT) - 0.172(a ) 2 + 
0.0106(AHA)(a) + 0.0697(AHA) 2 - 0.146(AOA)(a) - 0.04 16(AOA)( AHA) + 0.1 02 (AOA ) 2 - 
0.0694(OPTT)(a) - 0.00503(OPTT)(AHA) + 0.0151(OPTT)(AOA) + 0.01 73(OPTT) 2 (48) 

TTmax = 0.370 - 0.205(a) + 0.0307(AHA) + 0.108(AOA) + 1.019(OPTT) - 0.135(a) 2 + 
0.0141 (AHA)(a) + 0.0998(AHA) 2 + 0.208(AOA)(a) - 0.030 l(AOA)( AHA) - 0.226(AOA) 2 + 
0.353(OPTT)(a) - 0.0497(OPTT)(AOA) - 0.423(OPTT) 2 + 0.202(AHA)( a f - 

0.281 (AOA)(a ) 2 - 0.342(AHA) 2 (a) - 0.24 5 (AHA) 2 (AOA) + 0.281 (AOA ) 2 (AHA) - 
0. 184(OPTT) 2 (a) + 0.281(AHA)( a)(AOA) (49) 

X cc = 0.153 - 0.322(a) + 0.396(AHA) + 0.424(AOA) + 0.0226(OPTT) + 0.175(a) 2 + 
0.0185(AHA)(a) - 0.0701(AHA) 2 - 0.251(AOA)(a) + 0.179(AOA)(AHA) + 0.0150(AOA) 2 + 
0.0134(OPTT)(a) + 0.0296(OPTT)(AHA) + 0.0752(OPTT)(AOA) + 0.0192(OPTT) 2 (50) 

8. 4. Global Sensitivity Analysis 

While it is well established that variations in injector geometry can have a significant 
impact on performance and environmental objectives such as combustion chamber wall and 
injector face temperatures and heat fluxes [107], these sensitivities are seldom quantified. As 
shown in Vaidyanathan et al. [10], once surrogates become available, these sensitivities can be 


45 



computed using Sobol’s method (Section VI). 

In the context of the case study, a design variable is considered essential if it is responsible 
for at least 5% of the objective variability. Figure 13 shows the percentage of main factor (Si) 
contribution of different design variables to individual objectives. The variability of TF max is 
largely influenced by AHA and moderately by a (Figure 13a). The effect of the other design 
variables is marginal, suggesting that they are non-essential and in principle could be fixed. 
The variability of TW 4 is considerably influenced only by AHA (Figure 13b). TT max is 
influenced considerably by OPTT and marginally by a (Figure 13c). For X cc , AHA, AOA and 
a have considerable influence (Figure 13d). From the comparison of total sensitivity indices 
(S‘ otal ) with tiie main factors (Si) it is concluded that the contributions of the cross-interactions 
among the design variables to the objectives variability are negligible. 

The sensitivity information cited above can be used for screening purposes and therefore, 
ease the search for optimum designs. In fact, surrogate models that included only essential 
variables (with the non essential fixed at their mean values) were constructed and differences 
in the predictions were found insignificant when compared with those provided by the original 
(all design variables included) polynomial approximations. However, since the number of 
design variables in the case study is small, the original polynomial approximations were 
adopted. Note that the sensitivity analysis provides an insight into the physics of a design 
problem by highlighting the design features that govern the individual objectives. 

8.5. Surrogate-based Optimization 

Considering the problem of interest is one of multi-objective optimization, the optimal 
solutions are known as Pareto-optimal and their objective function space representation is 
called the Pareto-optimal front (POF). A solution is called Pareto-optimal (or efficient) 
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solution, if there is no other solution for which at least one objective has a better value while 
values of the remaining objectives are the same or better. Since surrogates for each of the 
objectives are available, any standard multi-objective optimization methods could be adopted. 
In the context of the case study, though, the results obtained by Goel et al. [11] using an 
evolutionary multi-objective optimization procedure (multi-objective genetic algorithm with 
e-constraint strategy) are presented. 

Extreme values of the Pareto-optimal front, that is, the results of optimizing (here, 
minimizing) each objective separately are shown in Table 5. A number one in parentheses 
indicates the objective selected for optimization. As expected based on physical grounds, the 
results show there is a strong correlation between TF max and TW 4 ; a correlation analysis shows 
a statistically significant correlation coefficient close to 1. As a result, the multi-objective 
problem can be addressed using only three of the original objectives ( TF max , TT max , X cc ). In 
addition, note that the geometrical characteristics of the optimal designs are aligned with the 
objectives, and as a result, the design that minimizes TF max corresponds to a shear-coaxial 
injector, that is, a = 0, and hydrogen flow area, AHA, and oxidizer post tip thickness, OPTT, 
at their maximum values; other objectives exhibit relative high values as expected for this type 
of injector. The results also show a correlation between TT max and X cc , although not as strong 
as the one for TF max and TW 4 . When TT max is minimized, X cc exhibits a low value, and vice 
versa. The optimal designs associated with the minimization of TT max and X cc correspond to 
impinging-like injectors with the hydrogen flow angle, a, at or near its maximum value and 
AHA and OPTT at their minimum values, yielding significantly shorter combustion lengths 
than those corresponding to shear coaxial designs. Both minimizations (TT max and X cc ) result 
in very high values for TF max . Considering that the optimization of each objective in isolation 
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resulted in alternative designs, and that objectives are in clear conflict (TF max , TT max \ TF max , 
X cc ) a full examination of the Pareto-optimal set is required. Figure 14 displays the two 
hundred and forty nine (249) Pareto-optimal solutions found. In order to analyze the POF a 
representative set of solutions (9) was identified using a hierarchical clustering algorithm. 
Values of the design variables and objectives for these designs are shown in Table 6. 
Graphically, these representative solutions are highlighted on the POF in Figure 14; note that 
the solutions are uniformly distributed along the POF. Box plots for the design variables and 
objectives in selected clusters (1, 3, 6, 9) are shown in Figure 15 and Figure 16, respectively. 
These plots highlight the variability of the design variables and objectives in each cluster. 

For cluster 1 it is seen that the value of a is fixed at 0 (shear coaxial injector) (Figure 15a) 
and AHA is fixed at 1 (Figure 15b). This suggests that in this cluster, the designs are sensitive 
to a and AHA, both of which reach their extreme values. The remaining two design variables 
AO A (Figure 15c) and OPTT (Figure 15d), vary over a range. It is observed that TF max is 
minimized (Figure 16a) where as X cc and TT max lie near their maximum and have little 
variability. Hence the designs in cluster 1 tend to minimize TF max and represent shear coaxial 
injector designs. The design variables AO A and OPTT do not influence TF max but affect the 
remaining objectives, X cc and TT max . Partial correlation coefficients are estimated to obtain the 
relationship between these design variables and objectives (Table 7). It is noticed that as AO A 
increases, X cc increases (Rcorr = 1.00) and TT max decreases (Ron- = -0.638). As OPTT decreases 
both X cc and TTmax decrease (Rcon- =1.00 for both). 

Similar observations can be made for clusters 3, 6 and 9. Figure 16a-c show that in 
different clusters, X cc and TT max decrease with increments in TF max (based on the median of the 
box plots). This highlights the trade-off between the objectives. Cluster 9 provides 
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information about the opposing trend as to what was observed in cluster 1. As objectives X cc 
and TT max are minimized (Figure 16b, c) an impinging injector design is obtained ( a ~ 1, 
Figure 15a) with TF max exhibiting high values (Figure 16a). The AHA is near minimum 
(Figure 15b) contrary to the design in cluster 1 where high AHA minimized TF max . The AO A 
has considerable variation which suggests that the variability of objectives, X cc and TT max are 
not largely affected by this design variable. Figure 15d shows that X cc and TT max are 
minimized for the minimum value of OPTT. Table 7 gives the partial correlation coefficients 
for the set of design variables and the objectives in each cluster which shows considerable 
variation. The partial correlation coefficients for the combinations left out are effectively zero. 

Finally, the Pareto fronts of the objectives known to be in conflict ( TF max , X cc ) and 
(TF ma x, TTmax) are shown. Figure 17a shows the relation between TF max and X cc . The O 2 post 
tip temperature, TT max is ignored. It can be seen that the POF in this case is linear over a large 
region. A small increase in the value of TF max (~ 10%) reduces the combustion length, X cc , by 
nearly 50%. Figure 17b shows the relation between TF max and TT max . The combustion length, 
X cc is ignored. It is obvious that the POF is non-convex. It is also seen that a small drop in the 
value of the face temperature (* 10%) can reduce the tip temperature TT max by nearly 60%. 
Hence at a small cost of TF max both X cc and TT max can improve considerably. 

9. Summary and Conclusions 

The fundamental issues that arise in the surrogate-based analysis and optimization (SBAO) 
of computationally expensive models such as those found in aerospace systems were 
discussed. The issues included the selection of the loss function and regularization criterion to 
identify the surrogates, design of experiments, surrogate selection and identification, 
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sensitivity analysis, and how to incorporate surrogates in optimization efforts when constraints 
are present. Some of these issues were demonstrated through the multi-objective optimal 
design of a liquid rocket injector. 

In the context of the case study, the SBAO approach demonstrated its effectiveness by 
solving a model for a critical problem in space propulsion: the multi-objective optimization of 
liquid injector rocket injector designs. 

New, potentially more effective methods and techniques are increasingly available to 
address some of the above-referenced issues, including: improvements on the commonly used 
orthogonal arrays and Latin hypercubes for design of experiments, merit functions for 
sampling that account for model appraisal information, alternative loss functions for modeling 
(e.g., e-sensitive) which has led to the development of more robust approximation schemes 
(e.g., support vector regression), the use of multiple surrogates to address the issue of model 
uncertainty, more efficient strategies for model selection and evaluation (e.g., bootstrapping), 
efficient global sensitivity methods for screening and general sensitivity evaluations, and 
multiple-surrogate based optimization strategies. 

SBAO is an active area of research, and has made substantial progress in addressing the 
analysis and optimization of a variety of aerospace systems. It has the potential of playing a 
vital role in the successful full-scale development of modem aerospace systems while 
effectively considering the competing needs of improving performance, reducing costs, and 
enhancing safety. 
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KERNEL FUNCTION 

ESTIMATION SCHEME 

AT(x; x t ) = (1 + x.x t ) d 

Polynomial of degree d (PRD) 

K{x\x i ) = \x-x l 

Linear splines (LSP) 

2 

K (jc; jc. ) = exp(- x-x t / ) 

Gaussian RBF (GRF) 


Table 1 Examples of kernel functions and related estimation schemes. 


METHOD 

SCOPE 

NATURE 

MODEL 

INDEPENDENT* 

Local Sensitivity 
Analysis 

LOCAL 

QUANTITATIVE 

YES 

Morris Method 

GLOBAL 

QUALITATIVE 

YES 

Scatter Plots 

GLOBAL 

QUALITATIVE 

NO 

Correlation 

Coefficients 

GLOBAL 

QUANTITATIVE 

NO* 

Variance-based 
Parametric (e.g., 
IFFD) 

GLOBAL 

! 

QUANTITATIVE 

NO 

Variance-based 

Non-Parametric 

GLOBAL 

QUANTITATIVE 

YES 


Table 2 A detailed review of different techniques in sensitivity analysis, and application 
examples, can be found in the works of Iman and Helton [68], Kleijnen [69], and Frey 
and Patil [70], and references therein. 

♦Independent from assumptions about the model being linear, etc. 

•[There are exceptions e.g., the Spearman rank-order correlation coefficient. 


t _ i_i ~ 

Minimum 


Maximum 


Variable - 

Actual 

Scaled 

Actual 

Scaled 

a 

0° 

0 

a 0 

1 

AHA 

Baseline 

0 

Baseline+25% 

1 

AOA 

Baseline-40% 

0 

Baseline 

1 

OPTT 

Yin 

0 

2x in 

1 


Table 3 Range of design variables (a is an acute angle in degrees and X is the thickness of 
OPTT in inches). 
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Case 

a 

AHA 

AOA 

OPTT 



TT max 


X 

1.0 


0.0 

0.0 

0.998 



-0.004 

Y 

0.0 

0.5 

0.5 

1.0 

0.285 

0.395 

0.923 

0.567 


Table 4 Independent and dependent variable (Objectives) for Cases X and Y from CFD 
computations (Normalized values shown). 



a 

AHA 

AOA 

OPTT 

TF max 



x cc 

1 

0 

1 

0.592 

1 


HS3B8I 

BIBBS! 


CFD 





-0.00207 

0.0656 


0.758 






0.110 

0.080 


0.250 

2 

0 

1 

0 

1 

0.0309 (0) 




CFD 





0.0910 

0.0461 

0.911 

0.568 

EHSR3I 





0.900 

0.570 

2.84 

2.93 

3 

1 

0 

1 

0 

BBHBKiB 











0.103 







0.0100 




4 

0.917 

0 

0 

0 

0.987 (0) 


0.182 (0) 


CFD 





0.981 

0.919 

0.119 

-0.004 






0.0800 

0.0800 

2.69 

0.120 


Table 5 Minimizing individual objectives. Value in parenthesis (1) indicates which objective 
function is minimized (Normalized values shown). 


Cluster 

a 

AHA 

AOA 

OPTT 

TF m ax 

Xcc 

TT max 

1 

0.000 

1.000 

0.842 

0.712 

0.0231 

1.090 

0.880 

2 

0.000 

1.000 

0.356 

0.587 

0.0276 

0.749 

0.890 

3 

0.000 

1.000 



0.0541 

0.750 

0.452 

4 

0.0939 

1.000 

0.000 

0.0146 

0.126 

0.453 

0.466 

5 

0.668 

1.000 

0.732 

0.000 

0.259 

0.681 

0.229 

6 

0.600 

0.670 

0.000 

0.000 

0.489 

0.264 

0.226 

7 

0.295 

0.108 

0.000 

0.354 

0.719 

0.129 

0.641 

8 

0.314 

0.0656 

0.000 

0.0554 

0.776 

0.0969 

0.357 

9 

1.000 

0.0140 

0.680 

0.000 

0.935 

0.138 

-0.0432 


Table 6 Objective function and design variables of nine (9) representative designs from the 
Pareto optimal solution. 
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Cluster 

Design Variable 

TF max 

Xcc 

TT max 

1 

AOA 

- 

1.000 

-0.638 

OPTT 

- 

1.000 

0.991 

3 

AOA 

- 

1.000 

- 

6 

a 

0.982 

-0.735 

-0.983 

AHA 

-0.999 

0.994 

-0.729 

9 

a 

0.877 

-0.203 

-0.769 

AHA 

-0.992 

0.983 

0.816 

AOA 

-0.977 

0.997 

-0.940 


Table 7 Partial correlation coefficients (Rcorr) of design variables vs. objectives for different 
clusters along the POF. 
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If necessary 



Figure 1 Key stages of the surrogate-based modeling approach. 



Figure 2 Anatomy of surrogate modeling: model estimation + model appraisal. The former 
provides an estimate of function / while the latter forecasts the associated error. 
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Figure 3 Multiple surrogates may be consistent with the data. For a given problem we oiay 
have preference for one or the other (e.g., additional information about the 
function / is available), but often which is the best surrogate is not clear a priori. 



Figure 4 A surrogate modeling scheme provides the expected value of the prediction 
E(f ) (solid line) and the uncertainty associated with that prediction illustrated 

here using a probability density function 9 . 



(a) Quadratic 


(b) Laplace 



(c) Huber (d) e-insensitive 


Figure 5 Alternative loss functions for the construction of surrogate models. 



Figure 6 A Latin Hypercube Design with N=6, K=2 for A" uniformly distributed on the unit 
square. 
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Figure 7 LHS designs with significant differences in terms of uniformity (Leary et al. 
[36]). 



Figure 8 The basic SMF algorithm. 
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Figure 9 Schematic of Hybrid Boeing Element (U. S. Patent 6253539). 




O z Flow Area 
(baseline-40%) 


H 2 Flow Area 
(baseline+25%) 


H 2 Flow Angle (0°-a°) 


0 2 Post Tip Thickness 
(x-2x in) 



Figure 10 Design variables and objectives of the single element rocket injector, (a) Design 
variables (b) Objective functions. 



adiabatic walk 



Figure 11 Simulation domain and boundary conditions of the injector flow model. 



Figure 12 Comparing the unrefined (336x81) {thicker lines} and refined (430x81) {thinner 
lines} grids. 
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Pareto optimal solutions $ 



Figure 14 Pareto optimal solution set and nine (9) representative solutions from the same 
set. The circles identify the representative of a specific cluster. 
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cluster cluster 

(c) (d) 

Figure 15 Box plots for the design variables in cluster 1, 3, 6 and 9. (a) a, (b) AHA, (c) 
AOA and (d) OPTT. 
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Figure 17 Pareto optimal front: (a) TF max vs. X cc (b) TF max vs. TT max . 
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With rapid progress made in employing computational techniques for various 
complex Navier-Stokes fluid flow problems, design optimization problems traditionally 
based on empirical formulations and experiments are now being addressed with the aid of 
computational fluid dynamics (CFD). To be able to carry out an effective CFD-based 
optimization study, it is essential that the uncertainty and appropriate confidence limits of 
the CFD solutions be quantified over the chosen design space. The present dissertation 
investigates the issues related to code verification, surrogate model-based optimization 
and sensitivity evaluation. 

For Navier-Stokes (NS) CFD code verification a least square extrapolation (LSE) 
method is assessed. This method projects numerically computed NS solutions from 
multiple, coarser base grids onto a finer grid and improves solution accuracy by 
minimizing the residual of the discretized NS equations over the projected grid. In this 
dissertation, the finite volume (FV) formulation is focused on. The interplay between the 
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concepts and the outcome of LSE, and the effects of solution gradients and singularities, 
nonlinear physics, and coupling of flow variables on the effectiveness of LSE are 
investigated. 

A CFD-based design optimization of a single element liquid rocket injector is 
conducted with surrogate models developed using response surface methodology (RSM) 
based on CFD solutions. The computational model consists of the NS equations, finite 
rate chemistry, and the k- s turbulence closure. With the aid of these surrogate models, 
sensitivity and trade-off analyses are carried out for the injector design whose geometry 
(hydrogen flow angle, hydrogen and oxygen flow areas and oxygen post tip thickness) is 
optimized to attain desirable goals in performance (combustion length) and 
life/survivability (the maximum temperatures on the oxidizer post tip and injector face 
and a combustion chamber wall temperature). A preliminary multi-objective optimization 
study is carried out using a geometric mean approach. Following this, sensitivity analyses 
with the aid of variance-based non-parametric approach and partial correlation 
coefficients are conducted using data available from surrogate models of the objectives 
and the multi-objective optima to identify the contribution of the design variables to the 
objective variability and to analyze the variability of the design variables and the 
objectives. 

In summary the present dissertation offers insight into an improved coarse to fine 
grid extrapolation technique for Navier-Stokes computations and also suggests tools for a 
designer to conduct design optimization study and related sensitivity analyses for a given 
design problem. 
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CHAPTER 1 

INTRODUCTION AND SCOPE 

Introduction 

With the advancements in computational technologies, different complex design 
problems can now be analyzed with reduced economical and computational cost as 
compared to an experiment. Additionally, computational models can be used as a 
substitute to study environments that are too hazardous to conduct physical experiments. 
These aspects have motivated the development of different tools in the areas of 
computational fluid dynamics (CFD) and multi-disciplinary optimization (MDO) which 
in turn promotes CFD-based design optimization studies. 

Computational fluid dynamics-based design studies can be broadly divided into 
two parts. The first part is the accurate modeling of the fluid flow problem and the second 
is the sensitivity estimation of the CFD solutions to the design variables and exploration 
of optimum designs. Solving fluid flow problems based on Navier-Stokes (NS) equations 
and other related models introduces challenges that include physical modeling 
uncertainty, geometric complexities, non-uniform and non-orthogonal meshing, disparate 
length scales and characteristics of the flow variables (such as velocity and pressure), and 
acceptable run time for engineering analysis and design. To obtain an accurate numerical 
solution after addressing all these issues is an elaborate topic of research. As pointed out 
in the review by Oberkampf and Trucano (2002), there are 5 sources of error in a CFD 
analysis, namely, insufficient spatial discretization convergence, insufficient temporal 
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discretization convergence, insufficient convergence of an iterative procedure, computer 
round-off and computer programming error. Much work has been done to estimate the 
first three errors. The remaining two are related to the field of computer programming 
and usually involve the verification of the code or the software. 

In order to ascertain the overall accuracy and identify uncertainties associated with 
physical modeling there is a strong need to develop techniques for accuracy assessment 
of a given numerical approach and mesh resolution. A priori and a posteriori error 
estimates are widely used to estimate errors and verify CFD solutions. A priori error 
estimation uses the information available before a solution is obtained based on the 
partial differential equation (PDE) that has to be solved and the initial and boundary 
conditions. In CFD analysis these estimates aid in the estimation of the stability and 
existence of a solution. On the other hand, an a posteriori error estimate uses the 
information of the computed solution. This allows the performance of a given numerical 
approach to be judged by estimating the errors arising out of the discretization of the 
PDEs. 

The issue of mesh resolution is of importance because a coarse grid is 
computationally economical but maybe inadequate in capturing the detail flow features 
whereas a fine grid, although it provides a well-resolved solution, increases the 
computational cost. Hence it is of interest to use a coarse grid and obtain an accurate 
solution. A well known method in this category is the Richardson extrapolation (RE). 
Details of this method can be found in Roache (1997) and Shyy et al. (2002). The 
Richardson extrapolation is based on eliminating the leading order error term based on 
the local truncation error analysis, implemented via solutions obtained on two or more 
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base grid levels. An attractive feature of this method is that it can be applied to any partial 
differential equation (PDE) without direct information about the solution procedure or the 
equation itself. On the other hand, the disadvantage of this method is that it assumes the 
order of the solution a priori, whereas generally the order of convergence of a CFD 
solution is space and solution profile dependent, and may not be consistent with that 
indicated by the local error analysis. As pointed out by Shyy et al. (2002), the flow field 
with high gradient in flow properties can deteriorate the performance of this method. 
Typically, RE is usually conducted on two base grids with one having twice the number 
of nodes as the other in each direction. For practical problems with complex geometries, 
it is often difficult to satisfy such a requirement. A very coarse grid may not efficiently 
capture the flow features for extrapolation and a very fine grid might increase the 
computational cost, thereby spoiling the whole purpose of extrapolation. 

Recently, Garbey and Shyy (2003) have proposed a least square extrapolation 
(LSE) technique. This method addresses some of the disadvantages noticed in RE. Unlike 
RE, LSE directly uses the discretized PDE model under consideration to estimate the 
degree of convergence. In this method residual is estimated on the fine grid onto which 
the base grid solutions are extrapolated based on the discretized equation and the specific 
computational scheme. Additionally the space dependency of the solution’s convergence 
order can be accounted for by assuming spatially dependent weight functions. Garbey 
and Shyy (2003) have shown that LSE is effective even when the base grids are not based 
on grid doubling or halving. Details of the LSE method will be discussed in chapter 2. 

In the area of MDO tools like design of experiment (DOE) (Myers and 
Montgomery, 1995) and response surface methodology (RSM) (Myers and Montgomery, 
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1995) have been shown to work well for a wide range of design problems. Hence it is of 
interest to use these tools to design devices that have complex geometries and flow 
features. Sensitivity analysis and optimization studies related to these devices require 
numerous function evaluations and the computational cost of CFD simulations restricts 
its use in such scenario. For such situation, RSM has been found to be very useful. A set 
of designs can be identified using DOE techniques and the solution for the designs 
obtained using CFD analyses. These solutions are then represented with response surface 
approximations (RS As) (polynomials or any other analytical function), which can then be 
used as the function evaluator. 

Sensitivity analysis helps the designer by providing the knowledge of the influence 
of design variables on the essential flow features. Global sensitivity analysis (Sobol, 
1993) is a method which has the potential for providing such information. The global 
sensitivity indices can be used to rank the design variables in order of their importance 
and also to identify the interaction between them. Based on these observations, 
nonessential variables can be fixed and the dimensionality of the design problem reduced. 

Most of the design studies require that more than one aspect of the solution 
(objective) be improved upon at the same time. There are different methods available for 
solving such multi-objective problems. Multi-objective optimization problems usually 
have numerous optimal solutions, known as Pareto optimal solutions (Miettinen, 1999). 
These solutions are such that no improvement is possible in any objective without 
sacrificing at least one of the remaining objectives. Hence these solutions are non- 
dominated. On the other hand if for a given solution there are other solutions where 
improvement in any objective is possible without sacrificing the remaining objectives, 
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that solution is said to be dominated. The function space of all the non-dominated 
solutions in the Pareto-optimal set is termed the Pareto-optimal front (POF). This set 
provides the designer with a clear idea of the trade-offs involved in a design study. 
Evolutionary algorithms (EAs) are popular global optimizers that have been used to find 
multiple Pareto optimal solutions in a single simulation run. Different multi-objective 
evolutionary algorithms (MOEA) have been proposed in literature (Deb, 2001). Insight 
into the trade-offs between different objectives can be obtained from the POF. 

In this dissertation different aspects of NS code verification and CFD-based design 
process are addressed and methods explored with the help of individual case studies. The 
LSE used for NS code verification is implemented on a lid-driven cavity flow with 
different boundary conditions and Reynolds numbers. The goal is to assess the 
effectiveness of the LSE technique for finite volume (FV) NS computations. 

The CFD-MDO coupled design study is addressed with the aid of a single element 
liquid rocket injector. The study is conducted using response surface methodology along 
with global sensitivity analyses and genetic algorithms. The overall goal is to estimate the 
global sensitivity and trade-offs involved in the design of a rocket injector. 

Scope 

The scope of the dissertation can be broadly divided into two, namely, 

• Investigation of NS CFD code verification. 

• CFD-based design optimization. 

The NS CFD code verification deals with the demonstration of LSE for numerical 
accuracy improvement and computational cost reduction with the aid of few widely used 
benchmark problems. The implementation of the method for complex flows involved in 
injector design is an issue of future study. The goal of the LSE related study is to 
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• Assess the effectiveness of the LSE technique for FV computations. 

• Evaluate the implication of solution gradients and singularities on the performance 
of LSE. 

• Address the issue of coupling between the pressure and velocity in the NS 
equations. 

The CFD-based design optimization study is addressed by modeling a single 
element liquid rocket injector. The goal is to integrate CFD with the optimization tools 
and obtain better design methodologies than those used in the past. While the ultimate 
goal is to analyze multi-element injectors, much of the detailed work in injector design 
can be done, or at least initiated, at the single element level. 

Design variables governing the life/survivability and performance of the injectors 
are identified, namely, H 2 flow angle, H 2 and 0 2 flow areas with fixed mass flow rates of 
fuel and oxidizer and 0 2 post tip thickness. The design objectives that are 
life/survivability indicators are the maximum temperature on the oxidizer post tip, the 
maximum temperature on the injector face and the combustion chamber wall temperature 
taken three inches from the injector face. The performance indicator is the length of the 
combustion zone. 

To facilitate the development of the present methodology, a baseline element 
design is needed. This baseline concept is generated by an empirical design methodology 
based on a specific set of propellant flow rates, mixture ratio and chamber pressure. The 
selected design variables are then varied based on this baseline design and the design 
space populated with the aid of a DOE technique. The prescribed CFD cases are executed 
and post processed to extract the required dependent variable data. This data are then 
used to generate a RSA for each objective in terms of the independent design variables. 
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Sensitivity and trade-off analyses for the design are then carried out using the data 
obtained from surrogate models of the design objectives, and the POF (multi-objective 
optima) generated by Goel et al. (2004) with the aid of a multi-objective genetic 
algorithm (MOGA) and a local search method (e -constraint strategy). The regions of the 
POF that represent different trade-offs among the objectives are obtained through a 
hierarchical clustering algorithm which will be discussed in chapter 3. This study is 
broadly divided into the following two parts. 

1. A sensitivity study is carried over the whole design space. The contribution 
of the design variables to the objectives variability is calculated using a 
variance-based non-parametric approach and correlations between objectives 
are investigated. 

2. Sensitivity analyses are conducted on clusters along the POF. Box plots are 
used to highlight the variability of the design variables and the objectives 
within a cluster. Additionally, the linear relationships between the design 
variables and the objectives are explored with the aid of partial correlation 
coefficients. 

The remaining chapters of this dissertation address each of the issues mentioned 
above in detail. In chapter 2 the different aspects of NS CFD code verification are given. 
In this chapter a literature survey is presented initially followed by the details of 
Richardson extrapolation (RE) and least square extrapolation (LSE) techniques. These 
techniques are implemented on a two-dimensional turning point problem and a laminar 
lid-driven cavity flow. The performances of the RE and LSE are evaluated and issues 
related to their implementation are identified. 

In chapter 3 different tools used for the design optimization study are presented in 
detail. Response surface methodology, DOE techniques, MOGA, optimization algorithm 
and sensitivity analyses tools are presented. Additionally a clustering method and a 
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visualization tool are presented which along with other tools mentioned above are used in 
chapter 4 to study the design problem in detail. 

In chapter 4 the past and present studies related to the design of a single element 
rocket injector are presented. The design variables and objectives of the current injector 
design are identified and optimization studies carried out for different design goals using 
surrogate models developed based on CFD solutions. Sensitivity analyses are done with 
the aid of surrogate models and sensitivity indices. The design trends are observed with 
the aid of results obtained from these studies and visualization techniques. 

In chapter 5 the dissertation is summarized and conclusions drawn. The unresolved 
issues are identified and scopes for future studies identified. 



CHAPTER 2 

NAVIER-STOKES CODE VERIFICATION 
In this chapter, the literature survey related to the work done on computational fluid 
dynamics (CFD) code verification is presented. Following this, details of Richardson 
extrapolation (RE) and least square extrapolation (LSE) techniques are presented along 
with case studies and results. 

Literature Review 

In the area of CFD code verification, there has been considerable research done to 
address the issues of grid sensitivity and quantification of uncertainties in the 
computations (Roache, 1997; Pelletier et al., 2001; Gu et al., 2001; Turgeon et al., 2001). 
The governing partial differential equations (PDEs) of fluid mechanics are either linear or 
nonlinear depending on the problem solved. The most frequently used computational 
methods for solving these PDEs are finite element (FE), finite difference (FD) and finite 
volume (FV) methods. To obtain an accurate solution using any of these discretization 
schemes, adequate grid resolution is essential since for a consistent scheme the computed 
solution approaches the exact continuous solution as the grid spacing tends to zero (Shyy, 
1994). The error arising out of inadequate grid distribution falls under the category of 
domain discretization error. At the same time the discretization operators used for the 
derivatives in the PDEs should be according to the order of variation of the solution over 
the domain. If an improper variation is assumed, then it gives rise to the equation 
discretization error. Hence the goal is to achieve a fine grid distribution over the domain 
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and a valid discretization of the derivatives in the PDEs without violating issues related to 
the stability of the computational model. Additionally, to obtain a solution for these PDEs 
valid boundary and initial conditions need to be specified which will define the flow 
field. 

The error estimation can be done before or after the CFD analysis is carried out. 

The former is called a priori error estimate and the later is referred to as a posteriori error 
estimate. In this dissertation the aim is to implement the a posteriori error estimates to a 
Navier-Stokes (NS) CFD code. Hence a literature survey of the available a posteriori 
error estimates is presented. As pointed out by Oden (2001), there are two broad classes 
of error estimation methods, namely, residual methods and recovery methods. In the 
residual method, the residual estimates the degree with which the approximate solution 
fails to satisfy the equations of the original problem (Babuska and Rheinbolt, 1978; 
Demkowicz et al., 1984). The residual estimates are used to define error norms which are 
largely used for controlling adaptive meshing procedures. The other type, referred to as 
the recovery method, attempts to enhance the computed solution during the post- 
processing step, like Richardson extrapolation (Roache, 1997; Shyy et ah, 2002) and least 
square extrapolation (Garbey and Shyy, 2003). 

There is a vast amount of literature available with respect to a posteriori error 
estimators for finite element approximation. A recent book by Ainsworth and Oden 
(2000) provides a survey of the available literature on such error estimators. Babuska and 
Rheinbolt (1978), Demkowicz et ah (1984), Zienkiewicz and Zhu (1987) are some of the 
earliest works in this area. These works concentrate on error estimation in finite element 
problems. Zienkiewicz and Zhu (1992a, 1992b) in a two-part paper introduced a recovery 
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technique which was an improvement on their previous effort. They achieved higher 
order accuracy by using this approach and also found it to be cost effective. Their 
approach uses a single and continuous polynomial expansion for the nodal values of the 
solutions (displacements, stresses) over a patch of elements surrounding the node of 
interest. A least square method or an L 2 projection is used to obtain this polynomial. The 
error in the computed solution is then measured by comparing it with the recovered 
solution using different forms of error norms. Ainsworth and Oden (1992, 1993a, 1993b) 
have presented an element residual method for finite element computations. This method 
has later been implemented by Jasak and Gosman (2001) for finite volume approach. 
Oden et al. (1994) extended the application of element residual methods to NS equations 
using a finite element approach. They designed error norms to measure the error in 
velocity and pressure for incompressible flow problems. 

Use of a posteriori estimator in finite volume approach has been of interest only 
over the past few years. Ilinca et al. (2000) have compared three different error estimators 
for finite volume solutions. They have compared a technique using Richardson 
extrapolation, a technique based on the difference of the computed solution, and a higher- 
order reconstruction obtained using a least square method and a technique which solves 
an error equation. Richardson extrapolation (Roache, 1994) uses the solutions on 
different grids with different levels of refinement, one finer than the other. Using the 
Taylor series representation of the discrete solution and combining the available multiple 
solutions, a higher order approximation of the desired variable is obtained. The drawback 
is to generate an adequately resolved solution on multiple grids for complex problems 
and the a priori assumption of the solution accuracy order. Solution reconstruction 
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(Barth, 1993; Ollivier-gooch, 1996) based on an approximation of the derivatives in the 
Taylor series expansion using a weighted least square method improves solution 
accuracy. This is then used to compare with the initial solution to obtain an error 
measure. The third method proposed by Zhang et al. (1997) estimates error using an error 
equation. In this method the error equation is derived directly from the governing 
equation and has a source term, which needs to be computed. This source term is 
evaluated using the flux differences at the cell interface. These methods were tested on 
subsonic, transonic and supersonic flows. Richardson extrapolation and the error equation 
method are shown to perform reasonably well as compared to the solution reconstruction 
method. 

Richardson extrapolation is based on eliminating the leading order error term in the 
assumed Taylor series expansion of the solution. This is done by combining solutions 
obtained on two different grids with (uniform) discrete spacings. The attractive feature 
about this method is that it can be implemented on a solution for any PDE without 
bothering about the solution procedure or the equation itself. On the other hand, the 
disadvantage of this method is based on the fact that it assumes a fixed order of accuracy 
(when using two grids) of the solution all over the computational domain. In practical 
fluid problems, due to numerical schemes and fluid physics involved, the order of 
accuracy is not uniform over the computational domain. Hence, the local truncation error 
based on Taylor series expansion may not represent the global accuracy of the solution. 
To address this, three grids of different resolution are required to estimate the order of 
accuracy of the solution. Additionally, problems involving turbulent flow features will 
result in sharp gradients in flow parameters in different regions of the domain. It has been 



13 


noticed by Shyy et al. (2002) that flow field with high gradient in flow properties can 
deteriorate the performance of this method. Although the most common use of RE is with 
two grids, one twice as fine as the other, this is not always essential. But for coarse grids, 
the assumption of monotonic truncation error convergence may not be valid thereby 
requiring fine base grids. For practical problems with complex geometries, it is tough to 
obtain two grids that satisfy such requirements. A very coarse grid may not efficiently 
capture the flow features for extrapolation, and a very fine grid might increase the 
computational cost, thereby spoiling the whole purpose of extrapolating. Another 
important disadvantage of RE is that the extrapolation does not preserve the 
conservativeness of the flow properties. 

The book by Roache (1998) has more information on different approaches to 
address some of these issues. But in its simplest form, RE is based on an a priori 
assumed order of convergence of the continuous solution. In the work of Garbey and 
Shyy (2003), basic properties of RE and its computational implications are presented in 
detail. They have presented complementary views that are asymptotic expansion for 
continuous function in a normed vector space and numerical approximation for discrete 
functions defined in a mesh. Comparing these views they have pointed out that it is 
difficult to achieve asymptotic order of convergence unless the numerical perturbations 
(arising out of imperfect convergence of the fluid solver) in the computations are small. 
Unless a grid is very fine, it is hard to reduce the perturbation error. Based on their study 
on the convergence order approximation and RE in relation to CFD problems, they have 
found that convergence order is space and grid dependent. Hence, they have concluded 
that there is no way to justify the efficiency of RE which uses a constant weight 
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formulation. Additionally, RE is easy to implement only when the grid spacing is 
uniform and the subsequent grids used in the extrapolation are obtained by grid doubling 
or halving. 

Scope 

The aim of this work is to investigate the performance of LSE method in NS 
computations (Vaidyanathan et al., 2004a). In particular, FV formulation is considered, 
which is different from the FD approach adopted by Garbey and Shyy (2003). The lid- 
driven cavity flow with different boundary conditions and Reynolds numbers (Re) is 
adopted as the main test problem. The goal of this study is to address the following 
issues: 

• To assess the effectiveness of the LSE technique for FV computations. 

• To evaluate the implication of solution gradients and singularities on the 
performance of LSE. 

• To address the issue of coupling between the pressure and velocity in the NS 
equations. 

It is noted that the geometric complexity plays a major role in practical flow 
computations. This aspect can have substantial impact on the performance of any 
extrapolation techniques. This issue will be addressed in future studies. 

In the following sections, the salient features of both RE and LSE techniques will 
be first presented. Following this, a brief discussion of the flow solver and the algorithms 
for the implementation of the LSE in relation to the NS computations will be discussed. 
Finally, their implementation on test cases and a discussion based on the results obtained 
will be provided. 
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Approach 

In the following, RE is reviewed first to help motivate the main topic of interest, 
namely the least square extrapolation. 

Richardson Extrapolation (RE) 

In essence, RE is based on an a priori assumed order of accuracy of the continuous 
solution. Considering a flow-field output, say a velocity component or pressure, of a CFD 
simulation on a grid size, h, and assuming second-order convergence a priori , the Taylor 
series expansion can be written as 

U(x) = u(x;h) + a 2 h 2 +a 3 h 2 +... (2.1) 


where U(x) is the exact solution and u(x;h) is the numerical solution based on h. Similarly 
the solution on a grid size, h/2 can be written as 


U{x)-u 
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Solving for U(x), using Eqs. (2.1) and (2.2) and eliminating 0(h 2 ) term and 
neglecting the higher-order terms leads to 
1 


U(x)=- 


4u 


+ O (A 3 ) 


(2.3) 


This is a second-order RE which results in a gain of order one. Similarly if the 
solution is assumed to have first-order convergence a priori and 0(h) term is eliminated, 
it results in a first-order RE with a gain of order one. 


U(x) = 2u x;^-u(x;h) + o(h 2 } 


(2.4) 


When in a practical fluids problem, the degree of convergence is not uniform over 
the computational domain, the order of extrapolation used can either improve the order of 



16 


accuracy or not affect it at all. If Eq. (2.4) was based on assumptions made in Eqs. (2.1) 
and (2.2), there would be no gain in the degree of convergence. Equations (2.3) and (2.4) 
can be generalized for two grids with spacing hi and h 2 , resulting in the following 
respective extrapolation scheme. 

u(x)= !^MtM!L S. ( 2 . 5 ) 


hfi(x,h,}-h,u(x,h.) _ 

u(x) = * K — y - ( 2 - 6 ) 

\K~h 2 ) 

where u(x,hj) represents the interpolated values from two coarse grid solutions, which are 
not necessarily based on grid doubling or halving. In implementing the extrapolation 
techniques, it is required that the coarse grid solutions be interpolated to the locations of 
the fine grid on which the extrapolated solution is obtained. The order of interpolation 
has to be such that it does not deteriorate the order of the computed solution. 

Least Square Extrapolation (LSE) 

In this dissertation the least square extrapolation method is used which, most 
importantly, estimates the order of convergence as the solution to a least square problem 
instead of assuming it a priori. Additionally, it accounts for the dependence of the order 
of accuracy on spatial coordinates by using spatially dependent extrapolating weights. 
The details of this method are presented below. In this approach, two coarse grid 
solutions, ui and U 2 , not necessarily based on grid doubling or halving, are combined 
linearly using a weighting function, a, which can be spatially dependent. The 
extrapolated solution is given as 

U =au l +(l-a)u 2 (2.7) 
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where u/ and u 2 represent the coarse grid solutions interpolated onto a fine grid, M 0 . 

Following Garbey and Shyy (2003), a modified Fourier expansion has been used 
for the weighting function, a. The weighting function, a is given by 

M 

a (jc) = a 0 + a x cos(x;r) + sin {(j- 1 (2.8) 

with (Xj,j = reals. The weighting function, oris dependent on the spatial co- 

ordinate, x for a one-dimensional problem. For two dimensional problem it is a function 
of spatial co-ordinates, x and y such that 

a(x,y) = a(x)a(y) (2.9) 

This modified Fourier expansion is ideal for rectangular domains. Different function may 
have to be used for more complex domains. The coefficients, ay are estimated by solving 
a least squares problem. For a given linear PDE, say 

L[u] = f (2.10) 

with u e (£ a ,|| || a ) and / € [E b ,\\ || fc j , its numerical approximation can be written as 

H=/* 

with u h e (£*,|| and f h e (£*,|| ||J . Based on this, the least square problem can be 
formulated as 

P a : Find a e /1(0,1) c L x such that 

L h [au x +(l-a)u 2 ]-f h (2.12) 

is minimum in L 2 (M C ,), where M 0 represents the fine grid to which the coarse grid 
solutions are interpolated. For a one-dimensional problem using two Fourier modes, Eq. 
(2.12) that has to be minimized can be rewritten as 
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L h [a 0 «, + (1 - a 0 )u 2 ] + . L h [a, cos {xn)u x + (l - a, cos(x;r ))u 2 ] ■ - /* (2.13) 

This approach can be generalized for nonlinear PDE problems via a Newton-like loop. 
For a given nonlinear problem, say, 

W[«] = / (2-14) 

linearization results in 

/(«)[v] = g (2-15) 

The least square problem then becomes 

J h {a k u l + (l-a*)w 2 )[a* +l M, +(l -a* +1 )« 2 ]-g A (2.16) 

which is minimum in L 2 (M 0 ). The iteration is started from the initial condition of = 0, 
until ||q:* +1 -a k | is less than some tolerance value. This method requires that the initial 

guess should not be too far from the final solution for convergence; i.e., the solution ui 
should be close to the grid solution on M 0 . In situations where J(u) is not available, a 
nonlinear set in a is obtained, which can be solved using nonlinear solvers. This is an 
issue of further study. In this initial implementation, a linearization is adapted to simplify 
the problem in hand. 

A point to be noted is that LSE uses interpolated functions on a fine grid, and the 
differential order, say n, of the PDE that is being solved impacts the choice of the 
interpolating scheme. The scheme used should give a result that should be smooth 
enough to be n-times differentiable. In the present study, cubic spline is used for 
interpolation of the solutions. 

The LSE technique is based on the minimization of the residual. However, the real 
goal of LSE is to maximize the solution accuracy, or equivalently, minimize the solution 
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error on the extrapolated grid. However, minimization of the residual is not necessarily 
the same as the minimization of the solution error. This aspect is explained in the 
following subsection. 

Error versus Residual 

To see the above mentioned point, let us consider only the linear problem, 
assuming the following form: 

AU = b (2.17) 

where A is a symmetric coefficient matrix of size nxn and U is the n x 1 exact solution 
vector of the linear system. Assuming u to be a non-exact solution of size n x 1, we have 
the following relationship: 

Au = b + r (2.18) 

where r represents the n x 1 residual vector of the system of equations. Defining e as the 
difference between u and U and subtracting Eq. (2.18) from Eq. (2.17) gives: 

As = -r (2.19) 

Denoting the eigenvectors of matrix A by a„ i = 1, ..., n such that ||a,|| = 1.0 and the 
eigenvalues as i = 1, ..., n the following relationships can be obtained: 

r = ( 2 - 20 ) 

/=1 

where i= l, ...,n are linear weights associated with each eigenvector and 

e = -A'‘r (2.21) 

which can be rewritten as 

<s = E(«,Mk (2.22) 

i-i 


The residual, F, is then defined as 



The least square approach aims at minimizing the I 2 norm of the residual, F; i.e., 


( n 

min (F) = min \ ^ |/; | 


( n 


= mm 


Ehf 


(2.24) 


v 1=1 


Now the adequacy of the non-exact solution, u, can be measure by the L 2 norm 
error measured between the u and U, which can be written as: 


error = E 


n 

=2>, 


(2.25) 


1=1 


Minimizing the error, gives 


n 

min (error) = min ^ |(a,. / )| 

l /*=! 


2A 


(2.26) 


From Eqs. (2.24) and (2.26) it is clear that minimization of F is equivalent to 
minimization of L 2 error norm {error) only when A,-’s are equal. Therefore it is clear that 
properties of the matrix A define the properties of 2.,’s and the relationship between F and 
L 2 error norm {error). The implication of the residual minimization in accuracy control 
will be examined along with the test problems. 

In the context of the current study, the solution error measure is defined as follows 


e: = 




1/2 


j= 1 1 

where n indicates the source of the non-exact solution, m indicates the source of the 
solution with respect to which the error is measured; i.e., the reference grid, (u.j ) 


(2.27) 


represents the extrapolated or interpolated solution (non-exact solution) onto the fine 
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grid, M 0 , at the (i,j) th node, ) represents the reference grid solution, at the (i,j)‘ h node 


and N is the number of nodes of the fine grid (M 0 ). The reference grid can be either the 
fine grid ( M 0 ) solution or a fixed grid solution (M„) fine enough to be considered as the 
benchmark solution. 

Flow Solver 

The governing equations are solved using the in-house code STREAM which is 
based on the pressure-based algorithm and finite volume approach (Thakur et al., 2002). 
The equations are solved on a multiblock structured curvilinear grid. The convective 
terms in the momentum equations are discretized using the second-order upwind scheme 
(Shyy, 1994) and the diffusion terms are discretized using a central difference scheme. 
The continuity and momentum equations in Cartesian coordinates are presented below. 
Steady state computations are carried out. 

Continuity: 

T-(/»)+T-(/»)-° < 2 - 28 > 

OX ox 

Momentum: 


— (pu(j)) + — (pv<f>) = — 
dx y ' dy KH ’ dx 


f 


df 


d<f> 


^ dx) + dyy dy j 


-S D 


(2.29) 


where p is the density, u and v are the velocity components in Cartesian coordinates, p is 
the laminar viscosity and Sp represents the pressure gradient term. In Eq. (2.29), <j> 
represents u or v component of the velocity. 


Figure 2.1 shows the collocated grid system. All the flow characteristics, namely 
pressure, velocity, density and viscosity are stored at the cell center (P). Carrying out the 
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volume integral on Eq. (2.29) and using the notations shown in Figure 2.1 for the nodal 
references, the following discretized equation is obtained for the momentum equations. 
Gptfip — &E$E a J.W — Sp + S c (2.30) 

where the a ’s represent the convective and diffusive fluxes at the respective grid 
locations, Sp is the source term arising out of the pressure gradient term and Sc is the 
source term arising out of higher order derivatives in the convective fluxes. More details 
of this equation can be found in Shyy (1994) and Shyy et al. (1997). The residual in the 
context of LSE is defined by subtracting the LHS from the RHS of Eq. (2.30). The flow 
solver is based on the SIMPLEC (Van Doormal and Raithby, 1984) algorithm where the 
pressure is updated by solving a pressure correction equation formulated out of 
combining momentum and continuity equations. For the current study, uniformly spaced 
Cartesian mesh is used since the case studies involve a square domain. 
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Figure 2.1 : Collocated grid and notation for 2D grid. 

In the implementation of LSE, a distinct equation is needed for each of the flow 
parameters that have to be extrapolated. Hence for pressure, the pressure equation is used 
which is obtained by substituting the velocity components, obtained from the momentum 
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equations, into the discretized (in the FV sense) continuity equation. The obtained 
pressure equation is given as 

A pPp ~ A E p E + A w p w + A N p N + A s p s + S (2.3 1) 

where the S is a function of the velocity components and the A’s are obtained by 
modifying the momentum and continuity equations (note that A ’s are different than the 
a’s in the momentum equation). The details of the pressure equation can be found in 
Patankar (1980). In summary, there are three equations to estimate the residuals required 
for LSE, namely Eq. (2.30) for u and v-component of the velocity and Eq. (2.3 1) for the 
pressure. The residuals in the context of LSE are defined by subtracting the LHS from the 
RHS of Eqs. (2.30) and (2.31). 

Algorithm for LSE 

As already mentioned there are issues that need to be addressed while 
implementing the LSE technique, namely, the non-linearity in the momentum equations 
(a’s have velocity components in them which make the equation non-linear) and the 
coupling of pressure-velocity so as to preserve the conservativeness of the flow properties 
in each cell. The pressure equation is linear in nature and therefore is simpler to work on 
provided that the source term is adequately represented. To proceed systematically, the 
Navier-Stokes problem is tackled in two steps. In the first step, only the pressure is 
extrapolated and in the second step the pressure-velocity coupled extrapolation is carried 
out. The algorithms are presented in the following sections. 

LSE for pressure equation only 

This step is to test the efficiency of the approach in extrapolating a given solution 
accurately; therefore the velocities obtained from the CFD solution on the finer of the two 
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coarse grids used for extrapolation are used for the source term in Eq. (2.3 1). The 
extrapolated pressure, plse, is defined as 

Plse ~ a Pi "*■ (l — ® ) Pi (2.32) 

where p { and p 2 are the interpolated pressure from the coarse grids to the fine grid using 
spline. The residual for the LSE method is then obtained by using the pressure equation, 
(Eq. (2.31)). For the preliminary analysis of the Navier-Stokes problem, only one mode 
(Eq. (2.8)) is used to obtain the pressure field. 

The algorithm for the LSE of pressure can be stated as 

• Define p LSE =a p ] + (1 -a)p 2 . 

• Input plse, U 2 > v 2 into Eq. (2.31). 

• Find a that minimizes the residual of Eq. (2.3 1). 

• Calculate plse- 

LSE for coupled pressure-velocity 

The momentum equation (Eq. (2.30)) is nonlinear in nature as the coefficients 
contain the velocity components as well. Hence, to linearize the system, a Picard-like 
iterative scheme (Ferziger and Peric, 1999) is adapted. The equation used for the velocity 
extrapolations comes from Eq. (2.30). 


a n d> n + l =a n f +l +a"f + ' +a n d> n+l +a n f +l -S P + S 

P ' f*LSE E E LSE W *LSE N T N LSE S T Slse ^ 


»+l 

C 


(2.33) 


where Sc +1 is defined based on the new velocity components. 

Now the velocity components and pressure are extrapolated sequentially. Again, 
only one mode (Eq. (2.8)) is used for the extrapolation of all the three parameters. The 
algorithm of the approach is given below 


START 
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• Define plse =ap l + (1 -a) p 2 . 

• Input plse, & 2 , v 2 into Eq. (2.3 1) {ulse and vlse are used from the 2 nd loop). 

• Find a that minimizes the residual of Eq. (2.31). 

• Calculate plse- 

• Define u L se =/ui + (1~Y) &2- 

• Input plse, ulse, v 2 into Eq. (2.33) with (f> replaced by u ( v L se is used from the 2 nd 
loop). 

• Find ythat minimizes the residual. 

• Calculate ulse- 

• Define vlse =P v, + (1-P)v 2 . 

• Input plse, ulse, v L se into Eq. (2.33) with <f> replaced by v. 

• Find p that minimizes the residual of the equation. 

• Calculate vlse - 

• Go back to START for the next iteration. Continue the loop imtil a, y and P 
converge within the specified tolerances. 

These algorithms are so designed that a direct access to the coding of the CFD 
software is necessary. The governing equations have to be modified to accommodate the 
estimation of the weights. 

Results and Discussion 

The least square extrapolation method is tested by implementing it on the following 

cases: 

• Linear, two-dimensional turning point problem (Figure 2.2) and 

• Laminar, lid-driven cavity flow (Figure 2.5). 

1 . variable lid- velocity 

2. constant lid- velocity 
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Two-Dimensional Turning Point Problem 

The linear equation that is solved is given below. 


&14 2 

eAu+a(x,y ) — = 0,xe(0,7r) 
dx 


(2.34) 


where £=0.1 and 


a(x,y) = x- 


n 


+ 0.3 


( 

y - -t 

V 2 y 


(2.35) 


H=0 


u=y(n-y) 

n 


u=-y(7T-y) 


Figure 2.2: 2D turning point problem (dimensions are n x x). 

The boundary conditions of the problem are defined such thatxj; e [0,7i]; u =y(n- 
y) at x = 0 and u = -y(n-y) at x = 7i; u = 0 aty = 0 and n (Figure 2.2). There is no velocity 
component along the y direction. This case has been tested by Garbey and Shyy (2003) 
where they have used a finite difference approach. In this work a finite volume approach 
is used. The problem is defined so that the transition layer of thickness e (where the 
velocity changes direction) centered on the curve a(x,y) = 0 is not parallel to x or y axis 
thereby making the problem two-dimensional. A sample solution is shown in Figure 2.3. 
In the finite volume implementation, central difference is used for the diffusion term and 
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first-order upwind is used for the convection term. Two modest base grids Mi and M 2 of 
sizes 23 x 23 and 29 x 29, respectively are chosen such that there is at least one or two 
grid points in the transition layer for the solution on grid M 2 . Both first-order and second- 
order RE are evaluated and compared with the performance of LSE. For this case study 
LSE is implemented with 4 Fourier modes for a in each spatial variable (instead of one 
mode as in the case of Navier-Stokes problem). The extrapolation is done onto fine grids, 
M 0 , ranging from size 41 x 41 to 101 x 101 in steps of 10. In Figure 2.4, subscripts, RE1 
and RE2 refers to RE of first- and second-order, respectively. The errors are represented 
as E b a which denotes the error in solution obtained using a certain method (RE or LSE) or 
CFD simulation on a given grid (represented as a) as compared to the solution obtained 
through CFD simulation on grid denoted by b (e.g. is the RMS error in u estimated 

using all the points in the computational domain, of M 0 , as compared to the interpolated 
solution from a finer grid (M„) of spacing h/2). 



Figure 2.3: Solution ( u ) contours for the 2D turning point problem, grid = 81 x 81. 
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Figure 2.4: Application: turning point problem with s = 0.1, x-axis is the number of 
points in each direction for the fine grid on which extrapolation is done. 

From Figure 2.4 it is seen that RE2 gives better results than RE1 for grid with poor 
resolution. This is because the transition layer is under-resolved and the dominant error 
comes from the second order diffusion term. As the grid gets finer, performance of 
is better than RE2 as the first order convective error dominates. For some intermedi ate 
grid (51 x 51), there maybe a cancellation between the convection and diffusion errors 
which might result in large improvements for RE as seen in Figure 2.4. Richardson 
extrapolation fails as the grid M 0 gets finer than N = 61. In all the cases, LSE predict the 
fine grid solution with an error less than that of the fine grid, M 0 , which is an 
approximation of the exact solution for N as large as 101. With the modest grid sizes of 
23 x 23 and 29 x 29, there are only one or two grid points in the transition layer which 
has a thickness of s= 0.1. Still, the LSE gives a solution, which is more accurate th an the 
fine grid solution on which the extrapolation is done resulting in a gain of more than one 


order of accuracy. 
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This problem is relatively easy to solve and is used to confirm the potential of LSE 
for FV implementation. Use of 4 modes showed better performance of LSE than using 
one mode. Increasing the number of modes beyond four provided marginal improvement 
in the extrapolated solution. Least square extrapolation is found effective for FV 
computations. Additionally it shows that LSE works well with flows having transition 
layers. For this case the resolution of the coarse grids is adequate. There are still issues 
like nonlinear equations, coupling of flow parameters, sharp gradients and singularities in 
the flow that needs to be addressed. To do so, lid-driven cavity flow is used and the 
details are given in the following section. 

Laminar Lid-Driven Cavity Flow 

Garbey and Shyy (2003) have addressed the lid-driven cavity flow problem in a 
vorticity-streamfunction formulation with a finite difference implementation. In the 
current study the complete NS equations are solved in 2D for the lid-driven cavity flow 
and the pressure along with u and v-components of the velocity field are extrapolated. A 
laminar incompressible flow computation is carried out in a square cavity of dimensions 
1 x 1 for Reynolds numbers 5.33 and 1000. This problem addresses the motion of the 
fluid in a square container induced due to the lid-velocity um in the positive x-direction 
(Figure 2.5). The remaining walls of the container have the no slip condition. 

Two different scenarios are considered. In the first case, the singularity at the end 
points of the sliding lid of the unit square cavity is avoided by choosing a variable speed 
as follows: 

u nd =-16x 2 (l -xf 


(2.36) 
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where uud is the w-component of the velocity at the lid. The v-component at the lid is zero 
and x varies from 0 to 1. In the second case, the lid velocity is kept constant, uuj = 1.0. 
This results in singularities in w-velocity at the lid comers where there is an abrupt 
change of uud from one to zero. 

U=Uud 

v=0 


u=0 

v=0 


. u=0 

* * v=0 

Figure 2.5: Two dimensional lid-driven cavity flow (dimensions 1 x 1). 

The Reynolds number for the flow is defined based on the mean lid- velocity. The 
w-velocity contours are shown in Figures 2.6 and 2.7 for Reynolds numbers 5.33 and 
1000, respectively. Some difference is noticed in the w-velocity contours at the upper 
comer regions of the cavity for both the Reynolds numbers. The negative variable uud 
results in a w-velocity profile which looks like the mirror image of the positive constant 
uud for Reynolds number 1000. The constant distribution of velocity results in the 
presence of singularity at the comers of the lid for the pressure as noticed from Figures 
2.8B and 2.9B. The magnitude of pressure near the comers for the case with variable um 
is about 10 times more than that for the case with constant uud (Figures 2.8 and 2.9). The 
presence of singularity will be an issue while implementing the LSE technique. In the 
later sections, the effect will be highlighted by means of the obtained results. 


u=0 

v=0 

y > 




0.3 



A B 


Figure 2.7: w-velocity contour for grid 171 x 171 and Re = 1000. A) M/,y = -16x 2 (l-x) 2 . B) 

uud= 1 . 0 . 

The positive and negative values of pressure identifies whether the pressure is 
acting on the lid or away from it, respectively. 
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Figure 2.8: Pressure along the lid for grid 121 x 121 and Re = 5.33. A) uud = -16x 2 (l-xf. 
B) uud = 1.0. 




Figure 2.9: Pressure along the lid for grid 171 x 171 and Re = 1000. A) uud = -16x 2 (l-x) 2 . 
B) uud = 1.0. 


The results for the following different scenarios are presented below: 

• LSE for Pressure only (Re = 5.33 and 1000) 

1. uud is variable. 

2. uud is constant. 

• LSE for Pressure- Velocity (Re = 5.33 and 1000) 
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1. uud is variable. 

LSE for pressure only 

The results are presented in two parts. In part (a) the case with variable lid- velocity 
is presented and in part (b) the case with constant lid-velocity is presented. In both the 
parts the results for Re = 5.33 are presented first followed by the results for Re = 1000. 
For Re = 5.33, the two selected base grids which captures the flow features adequately 
are 61 x 61 and 71 x 71. The extrapolated pressure fields are compared to the solution 
obtained using a fixed fine grid (M„) of size 121 x 121. The extrapolated fine grids (M 0 ) 
vary between sizes 81x81 and 121x121. For Re = 1000, the two coarse grids are 1 1 1 x 
1 1 1 and 121 x 121 and the fixed fine grid (M n ) for comparison is 171 x 171. The 
extrapolated grids (M 0 ) are between 131 x 131 and 171 x 171. Only one Fourier mode 
(Eq. (2.8)) is used for the extrapolation of pressure. The results are compared with results 
obtained from Richardson extrapolation. 

Variable lid-velocity (Re = 5.33). In Figure 2.10A the comparison of extrapolated 
solutions on grid M 0 (81 x 81, ...) with the CFD solutions on the same grid (M 0 ) show 
that LSE approach performs better as expected. The error increases as the grid gets finer 
since resolution of the base grid solutions is not adequate for extrapolating onto very fine 
grids. From Figure 2.1 0B, the true error can be estimated where the extrapolated 
solutions on any grid M 0 is compared with the fixed fine grid ( M „ ) solution which 
represents the “correct” solution, obtained on a substantially finer grid, for the given 
Reynolds number. As expected it is seen that as the grid gets finer, LSE outperforms 
other methods and attains an asymptotic convergence limit. The improved performance 
of LSE for finer grid is expected as the minimization of the residual should satisfy the 
governing equations more accurately. 
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Figure 2.10: Comparison of interpolated and extrapolated pressure with CFD solutions 
for Re = 5.33 with variable lid- velocity. A) Extrapolated pressure compared 
with CFD solution on the extrapolated grid (M 0 ). B) Extrapolated pressure 
compared with CFD solution on fixed fine grid (M„), 121 x 121. (E” : n refers 
to the source of the solution which is compared with the reference solution 
identified by m). 

To follow up on the previous discussion, the variation of the L 2 norm of the 
residual, F, and the L 2 norm error of the extrapolated solution on grid M 0 as compared to 
the fixed grid 121 * 121, with respect to a is examined. Figure 2.1 1 shows that the trend 
of both the accuracy measures are similar but based on minimization of F (Figure 2. 1 1 A) 
LSE picks a value of -1.48 for abut based on L 2 norm error (Figure 2.1 IB) a value of - 
1.6 for a is a better choice. Although it does not identify similar ds that minimize both 
measures, they are close to each other. More assessment will be presented for the flow 
with a constant lid velocity. 

Variable lid-velocity (Re = 1000). The performance of LSE is compared with 
other methods for Re = 1000 and the obtained results are shown in Figure 2.12. The base 
grids used for the extrapolation are 111 x 111 and 121 x 121. It clearly shows in Figure 





35 


2.12A that for a given fine grid (M 0 ) LSE performs better. Similar to the low Reynolds 
number case, the LSE solution improves as the grid gets finer in terms of the true error 
measure (Figure 2.12B). 


a vs F (Mo= 101) 
ocre 2 = -2.77, log, 0 (F) = -9.945 



A 


a vs E m „(Mo= 101) 
ccre2 = -2.77, log(o(E Mn RE2) = -3.588 



B 


Figure 2. 1 1 : Least square L 2 norm error and error in extrapolated pressure on grid 101 x 
101 for different a for Re = 5.33 and variable lid- velocity. A) Least square L 2 
norm error (F). B) L 2 norm error ( E ” ) (comparing extrapolated pressure with 
CFD solution on grid 121 x 121 (M„)). 

This exercise with variable lid-velocity and two different Reynolds numbers gives 
us incite into few important aspect of LSE. 

• The resolution of the coarse grid has to be adequate for efficient extrapolation. 
Hence for higher Reynolds number, resolution of the base grids is higher. 

• LSE performs better for a considerable range of Reynolds numbers (5.33 and 
1000), thereby showing its efficiency for a wide range of laminar and 
incompressible flow regime. 

• LSE performs better as the grid gets finer since it satisfies the governing equation 
more accurately over the computational domain. 

• The minimization of the residual does not guarantee the minimization of the 
solution error. 





36 



Figure 2.12: Comparison of interpolated and extrapolated pressure with CFD soluti ons 
for Re = 1000 and variable lid-velocity. A) Extrapolated pressure comp arec j 
with CFD solution on the extrapolated grid (Mo)- B) Extrapolated pressu re 
compared with CFD solution on grid 171 x 171. (£" : n refers to the sou rce 0 f 
the solution which is compared with the reference solution identified by m \ 

Constant lid-velocity (Re = 5.33). At this point it would be interesting to take a 
look at a case where the lid- velocity is uniform to evaluate the efficiency of LSE when 
singularities are present in the solution domain. 

Initially the LSE is implemented on the complete solution domain. In Figure 2.13, 
the sensitivity of L 2 norm of the residual, F and the L 2 norm of the solution error, £* ^ 

with respect to a is shown. LSE identifies a = -0.48 as the optimum value, based on the 
minimization of F (Figure 2. 13 A). But as seen from Figure 2. 13B this value of a is f ar 
from the value required to obtain an accurate solution which is about -4.0. 

To investigate the impact of the singularities, the LSE is carried out over a reduced 
solution domain. Figure 2.14 shows the reduced domain over which the LSE is 
implemented. The extrapolation is carried over a domain with x = 0-> 1 and y = O-^q. 95 . 
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a vs F (Mo = 101) 
ctRE 2 = -2-77, log 10 (F) = -6.282 



aVsE m „ (Mo ■ 101) 

(xre 2 = -2.77, log 10 (E% K ) = -2.365 
aisE = -0.48, log 10 (E M " LSe ) = -2.054 
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Figure 2.13: Least square L 2 norm error and error in extrapolated pressure on grid 101 x 
101 for different a for Re = 5.33 and constant lid-velocity (full domain). A) 
Least square L 2 norm error (F). B) L 2 norm error ( E” ) (comparing 
extrapolated pressure with CFD solution on grid 121 x 121 ( [M n )) . 



Figure 2. 14: The shaded domain is used for LSE. A region between y = 0 0.95 for all * 

is used for extrapolation. 

As shown before in Figure 2.13 for the whole domain, the sensitivity of least 
square error norm, F and the L 2 norm error ( £” ) in comparison to the CFD solution with 
respect to a for the reduced domain, of grid 101 x 101, is shown in Figure 2.15. It is seen 






38 


that LSE identifies a = -0.48 as the optimum value for a, based on the minimization of F 
(Figure 2. 15 A). But as seen from Figure 2.15B this value of a is not as far from the 
optimum value of a based on the minimization of E ” . The suggested value of a is 

approximately -1.3 as compared to the previous approximate value of -4.0. This suggests 
that a sub-domain which is adequately far from the region of singularity improves the 
performance of LSE. A point to be noted is that during interpolation of data from one 
grid to the other grid, the comer singularity tends to get smeared onto the final grid. Some 
method of treating or accounting for this might have to be looked into. 



A B 


Figure 2.15: Least square L 2 norm error and error in extrapolated pressure on grid 101 x 
101 for different a for Re = 5.33 and constant lid-velocity (reduced domain). 
A) Least square L 2 norm error (F). B) L 2 norm error ( E ” ) (comparing 
extrapolated pressure with CFD solution on grid 121 x 121). 

Figure 2.16 shows the comparisons of different extrapolation techniques. The 
errors shown are based on the reduced domain. There are oscillations noticed in the 
convergence rates as the grid gets finer. This is due to the fact that the region for different 
grids need not always lie between y = 0 -> 0.95, since the nearest grid point may be 
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below or above. This results in different levels of improvement. The present domain 
reduction strategy suggests that when a domain farther away from the singularities is 
considered, the performance of LSE seems better. Of course, this is a very simplistic 
solution and may not necessarily be able to handle other flow problems with more 
complexities. In Figure 2. 16 A, the extrapolated solutions are compared to the fine grid 
(M 0 ) solutions, which show that LSE outperforms the others. In Figure 2.16B, the 
comparison of the extrapolated solution is made with a fixed fine grid solution of size 
121 x 121. The solution on this grid is considered as the benchmark solution. Clearly 
LSE performs well as compared to RE. 


E m n VS N 
Re = 5.33 



Grid Points (N) 

A 



B 


Figure 2.16: Comparison of interpolated and extrapolated pressure with CFD solutions 
for Re = 5.33 with constant lid- velocity. A) Extrapolated pressure compared 
with CFD solution on the extrapolated grid {M 0 ). B) Extrapolated pressure 
compared with CFD solution on grid 121 x 121 (M n ). (E":n refers to the 
source of the solution which is compared with the reference solution identified 
by m) 

Constant lid-velocity (Re = 1000). LSE is implemented in the same reduced 
domain as shown in Figure 2.14. The coarse grids used are 1 1 1 x ill and 121 x 121 and 
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the fixed fine grid (M„), which is the representative of the correct solution, is 171 x 171. 
Figure 2.17, shows the comparison of different extrapolation schemes. Similar trends are 
noticed for the reduced domain as before. 



Figure 2.17: Comparison of interpolated and extrapolated pressure with CFD solutions 
for Re = 1000 with constant lid- velocity. A) Extrapolated pressure compared 
with CFD solution on the extrapolated grid (M 0 ). B) Extrapolated pressure 
compared with CFD solution on grid 171 x 171. ( E ™ : n refers to the source of 
the solution which is compared with the reference solution identified by m) 

This exercise helps gain an idea of shortcomings related to LSE technique. It shows 
that when singularities are present, the minimization of the L 2 norm of the discretized 
governing equation residual, F, need not minimize the solution accuracy measure, E ” . It 
depends largely on the matrix of equations used for the least square implementation. 
When the singularity is removed and the LSE is carried out on the reduced domain, the 
minimization of F tends towards the minimization of E” and the performance of LSE is 


improved. 
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LSE for coupled pressure-velocity computations 

The coupled algorithm is implement on a lid-driven cavity flow with Re = 5.33 and 
1000 and the variable lid- velocity as given in Eq. (2.36) is used. The geometry is same as 
before (Figure 2.5). The coarse grids {Mi, Mi ) , fine grids (M 0 ) and the fixed grids (M„) 
are all the same as those used for pressure extrapolation. It is noticed that about 3-4 
iterations of this algorithm results in converged values of a, y and with about 3-4 inner 
iterations for the velocity components. Only the error compared to the finest grid (M„) is 
shown since the goal is to obtain a solution as accurate as the exact solution. For clarity, 
only the RE2 and LSE solutions are shown in Figures 2.18 and 2.19 since they performed 
better compared to the rest. 

Re = 5.33. The E” measure for pressure, w-velocity and v-velocity are shown in 

Figures 2. 18 A, 2.18B and 2.18C, respectively. Single mode (constant, ao ) is used for the 
extrapolation of all the flow parameters. The LSE performs well in extrapolating all the 3 
flow parameters onto the finest grid. 

Re = 1000. In Figure 2.19, the LSE solution is compared with the solution of 
second order RE. There is considerable difference in the trends as compared to the low 
Reynolds number case. But as the grid gets finer, the performance of LSE is better than 
RE, suggesting that LSE gives a solution closer to the exact solution than RE. 

This exercise gives the following information: 

• The coupled algorithm is effective in treating all flow variables: pressure, u- and v- 
velocity, adequately, and works well for a wide range of Reynolds number. 

• It is surprising to notice that for the coupled algorithm, the improvement noticed in 
LSE as compared to RE2 seems to be more for the higher Reynolds number case 
t han the lower. Additionally the trend of w-velocity is reversed for Re = 1000 and 
as the grid gets finer, the performance of LSE is poorer in terms of accuracy. This 
requires further study. 
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Figure 2.19: Extrapolated pressure and velocity compared with CFD solution on grid 171 
x 171 ( M „ ), Re = 1000 and variable lid-velocity. A) Pressure. B) w-velocity. 
C) v-velocity. 







CHAPTER 3 

GLOBAL SENSITIVITY AND OPTIMIZATION TECHNIQUES 
In this chapter detail of the surrogate modeling method namely response surface 
methodology (RSM), global sensitivity approach and multi-disciplinary optimization 
(MDO) tools like design of experiments (DOE) and genetic algorithm are presented. 
These tools are applied in the design of a single element injector of a liquid rocket engine 
and the obtained results are presented in chapter 4. 

Response Surface Methodology (RSM) 

The approach of RSM is to perform a series of experiments, based on numerical 
analyses or physical experiments, for a prescribed set of design points (using DOE), and 
to construct a global approximation of the measured quantity (objective) over the design 
space. Figure 3.1 shows the three steps involved for a single-objective optimization study 
which is dependent on two design variables. This is a general case and can be extended 
for multi-objective optimization study as explained later in this chapter. Step I (Figure 
3.1) identifies the designs and the corresponding values of the objective for each design 
over the selected design space. In step II (Figure 3.1) a response surface approximation 
(RSA) of the objective is generated, which is a function of the two design variables. 
Finally step III (Figure 3.1) involves the search for the optimum design using the 
generated RSA for function evaluations. In this dissertation regression analysis is used to 
construct polynomial-based RSAs of assumed order and unknown coefficients. The 
number of coefficients to be evaluated depends on the order of the polynomial and the 
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number of design parameters involved. For instance, a second-order polynomial of N 
design variables has (N+l)(N+2)/(2!) coefficients. A cubic model has 
(N+l)(N+2)(N+3)/(3!) coefficients. In this dissertation, RSAs are constructed using JMP 
(2002), a statistical analysis software that provides a variety of statistical features in an 
interactive format. 

In the practical application of RSM, it is necessary to develop an approximate 
model (RSA) for the true response or the objective. The second order (quadratic) RSA is 
used most frequently as it is the most economical one. Such an approximation for a 
response variable (objective) y with k regressors can be written as 

y = Po+Ya P> x i + Z P» x i + Z Z Pa x i x ) +£ ( 3 -0 

1=1 1*1 i'=l 7=2 

The above equation can be written in matrix notation as follows 
y = Xp + e (3.2) 

where y is the (nx 1) vector of observations, X: (nxn p ) matrix of the levels of the 
independent variables, fi : (n p x 1) vector of the regression coefficients, e: (nx 1) vector of 
random error, n: the number of observations, and n p : the number of terms in the RSA. 

The purpose is to find the vector of least square estimators b that minimizes 

Error = Jf, 2 = e T e = (y - X P) r (y - X p) (3.3) 

1=1 


which yields the least squares estimator of P 
b = (X T X)-'X T y 


(3.4) 
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step I: Define and populate the space 


step II: Generate RSA. 


step III: Search for the optimum 




Figure 3.1: Schematic of the procedure for global design optimization (Papila, 2001). 
The quality of different RSAs can be evaluated by comparing the adjusted rms- 


error a a (Myers and Montgomogery, 1995) value that is defined as: 



where e, is the error at / th point of the data used for fitting. 

The accuracy of the RSA in representing the design space is measured by 
comparing their predictions to the actual values of the objective at test design points, 
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different from those used to generate the fit. The prediction rms-error <7 for the test set is 
given by: 


(3.6) 

In this equation £, is the error at the z th test point and m is the number of test points. 
The test points can be selected using the cross-validation technique (Papila, 2001), which 
is an established method for estimating the prediction accuracy. It is usually performed 
using either a number of random test/train partitions of the data, or a &-fold cross- 
validation (Mullin and Sukthankar, 2000). In &-fold cross-validation, the data points are 
divided in k equally sized subsets. One of the subset is kept for testing and the remaining 
k - 1 are used to fit the RSA. This is done for all the k subsets and the average error of all 
the k RSAs is computed. The fitting £-1 sets which has the minimum prediction rms- 
error, a, is selected and the remaining subset used as the testing set. 

In certain cases, especially when higher order polynomials are used, the data 
available for RSA may not be enough to spare some for testing the fitted model. Hence, 
an alternate method to estimate the performance of the RSA is to compute the PRESS 
rms-error. This method was proposed by Allen (1971, 1974) and it computes a sum of 
squares of the residuals. The residual is obtained by fitting a RSA over the design space 
after dropping one design point from the fitting set and then comparing the predicted 
value of the RSA for that point with the expected value. This is done for every point in 
the design space. The PRESS rms-error is given by 



PRESS 


1 


Eb .-£] 2 


/=! 


n 


(3.7) 
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where y,- is the expected value, y ( . is the value predicted by the RSA for the I th point, 
which is excluded while generating the RSA and n is the number of design points. If this 
value is close to the adjusted rms-eror, cr a , then the model performs well. 

A measure of the variability in an objective accounted by its RSA is given by the 
coefficient of multiple determination (Myers and Montgomogery, 1995), R , given as 

r2=1 _^e_ ( 3 . 8 ) 

SS 

yy 

n 

where SSe is the sum of squares of the residuals (=^(y,. -y,.) 2 ) where y is the predicted 

»=i 

value by the RSA. SSyy is the total sum of squares about the mean given by 

ss„=ss s +ss R =Y,(y,-y? 

i=J 

where y is an overall average ofy,. SSr is the sum of squares due to regression 


n 

(= ^ ( j>. -y ) 2 ) where y is the overall average of y,-. Since R 2 increases as more terms are 

i=i 

added in the RSA, the adjusted coefficient of multiple determination R 2 , a normalized 
measure is usually preferred. It accounts for the degrees of freedom in the model and is 


given by 



SS E /(n-n p ) 
SS yy /(n- 1) 


= 1 - 



(3.10) 


For a good fit, the value of R] should be closer to one. 

The polynomial-based RSA is effective in representing the global characteristics of 
the design space. It can filter the noise associated with design data. The linearity of the 
polynomial-based RSA allows us to use statistical techniques known as design of 
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experiments (DOE) to find efficient fitting sets. On the other hand, depending on the 
order of polynomial employed and the shape of the actual response surface, the RSA can 
introduce a substantial error in certain region of the design space. 

Design of Experiments (DOE) 

It is well established that the predictive capability of RSM is greatly influenced by 
the distribution of sampling points in design space (Unal et al., 1997, 1998). Design of 
experiment is required to select designs for RSA that minimizes the effect of noise. There 
are different types of DOE techniques in the literature as reported by Haftka et al. (1998). 
One of the most conservative DOE techniques available in literature is the full factorial 
design (JMP, 2002). Unal et al. (1996, 1997) studied response surface modeling using 
orthogonal arrays (O A) in computer experiments for reusable launch vehicles and 
illustrated that using this technique can minimize design, development, test and 
evaluation cost. Unal and Dean (1995) studied a robust design method based on the 
Taguchi method (Unal and Dean, 1991; Dean, 1996) to determine the optimum 
configuration of design parameters for performance, quality and cost. They demonstrated 
that using such a robust design method for selection of design points is a systematic and 
efficient approach for determining the optimum configuration. The full factorial design 
and the OA are explained below. 

In full factorial design the range of each design variable is divided into one or more 
intervals, which mark the number of levels. These intervals are usually, evenly spaced. 
All the possible combinations of the levels of all the design variables give the design 
points in the design space. For example, for three variables with three levels each there 
are totally 27 design points (Figure 3.2). 
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Figure 3.2: Full factorial 3-level design (dark circles refer to points in the foreground and 
light circles refer to points in the background). 

An OA is a fractional factorial matrix that assures a balanced comparison of levels 
of any factor or interaction of factors. Consider A, a matrix with elements of aj where j 

denotes the row (j = 1,2... n r ) and i denotes the column (i = l,2...n c ) that aj belongs to, 
supposing that each aj e Q = {0, 1 . . .q- 1 } . If the columns representing each design 
variable are mutually orthogonal, matrix A is called an orthogonal array. It is of strength t 
<n c if in each n r -row-by-t-column sub-matrix, the q‘ possible distinct rows occur X (index 
of the array) times. Such an array is denoted by OA(n r ,n a q,t ) by Owen (1992). Since the 
points are not necessarily at the vertices, the OA can be more robust than the full factorial 
design in interior design space and is less likely to fail the analysis tool. Based on the 
DOE theory, OA can significantly reduce the number of experimental configurations 
(Papila, 2001). In this dissertation OA is used to obtain the design for the optimization 
study. 
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Sensitivity Analysis 

Unlike local sensitivity where the partial derivatives are used to locally estimate the 
sensitivity of an objective to a specific design variable, global sensitivity allows the study 
of overall model behavior. To understand the concept, assume a surrogate model (in this 
dissertation a RSA) of an objective, J{X) as a function of the design variables, X, where X 
is the vector of the design variables xi, X 2 , ..., x„, scaled between zero and one. This 
surrogate model (in this dissertation, polynomials of n* order) represented in Eq. (3.11) 
has to be square integrable. 


f{,X ) — fo + ^ fi (*, ) + ^ fij (-*, , JCy )+ ... + f\ 2 ...n (*1 ’ X 2 >"’> X n ) 


(3.11) 


«</ 


An analysis of variance (ANOVA) (Archer et al., 1997) representation of Eq. (3.1 1) 
can be obtained by imposing the following condition 

=0 (3.12) 

0 


for k = ii , ..., i s . This condition is essential for the uniqueness of Eq. (3.11). This results 
in a set of summands of different dimensions. Each of the components is responsible 


for the interactive contribution of design variables ,...,x fj to the variability of the 

objective flX) over the n-dimensional unit hypercube design space. 

Integrating Eq. (3.1 1) over all the design variables gives 
j/(x)n<& =/„ (3.13) 

Integrating Eq. (3.1 1) over all the design variables except x, gives 

f/«n d ** = f & + ^ ( x i ) < 3 - 14 ) 

k*i 


from which/(x,) is obtained. 
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Integrating Eq. (3.11) over all the design variables except jc, and xj gives 
|/Wn^=/o +fiM+fj(xj)+f,j(x lt x j) (3.15) 

k*i,j 

from which fjixuxj ) is obtained. This can similarly be extended to obtain the remaining 
summands offlX). The imposed condition ensures that the summands are orthogonal in 
nature and since/(X) is square integrable the summands are too. Therefore the partial 
variances can be calculated as 


D, i = f fj 2 1 dx ...dx i 

l i— l s J*' 1 !— 1 ! *1 l t 


and the total variances is 


(3.16) 


D= jf 2 (x)dx-f 0 2 (3.17) 

Therefore it can be shown (Sobol, 1993) that 


n n 

0 = L 2X..,. 


ij <...</, 


(3.18) 


Partial variance gives a measure of the spread of the function due to one of the 
independent variables. Total variance gives a measure of the spread of the dependent data 
due to all the independent variables. Sensitivity is a measure of the contribution of an 
independent variable to the total variance of the dependent data. Sobol (1993) had 
proposed a variance-based non-parametric approach to estimate the global sensitivity 
using Monte Carlo methods. To calculate the total sensitivity of any design variable, say 
xi, the design variable vector X is divided into two complementary subsets, xi and Z 
where Z is a vector containing^, X 3 , ..., x„. The purpose of using these subsets is to 
isolate the influence of x, on f(X) variability from the influence of the remaining design 
variables included in Z. The total sensitivity index for*/ is then defined as 
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(3.19) 

where 

D? = Z>,+Z>,„z (3.20) 



D x is the partial variance associated with x; and D x Z is the sum of all partial variances 

associated with any combination of the remaining variables representing the interactions 
between xj and Z. Similarly the partial variance associated with Z can be defined as Dz- 
Therefore the total objective variability can be written as 

D-D^D Z +D^ (3.21) 

In cases where the/(X) is an analytical function, the multidimensional integrals for 
computing the partial variances can be evaluated analytically. However, the Monte Carlo 
approach is applicable under general conditions (e.g., any model, design under 
uncertainty) and has been adopted in this dissertation. The designs for the Monte Carlo 
approach are selected without any preference of one design over the other from the 
uniformly distributed design space of unit sides. Hence the design variables are 
uncorrelated which is essential for the effective implementation of the method proposed 
by Sobol (1993). The variance estimates can be obtained using the following expressions: 

(322) 

where fo is the mean of the objective. 


* ’ y=l 


(3.23) 


JV j ={ 


(3.24) 



where the terms (xj, Z) and (xi ’, Z’) represent random sample designs. Equations. (3.23), 
(3.24) and (3.25) give the estimates for D Xi , D and Dz, respectively and D z can be 

calculated using Eq. (3.21). Once these estimates are known the total sensitivity index of 
the objective variability with respect to a given design variable can be obtained using Eq. 
(3.19). The influence of a design variable to an objective variability without accounting 
for any of its interactions with other variables is denoted as a main factor index and given 
as 

S„=DJD (3.26) 

Each pair of random samples requires three different objective function evaluations 
(e.g .,J{x], Z ),/ (xj, Z’) and / (xj Z)). The mean and the total variance of an objective 
need to be estimated only once, and only two Monte Carlo integrals per design variable 
are necessary to compute the main factor and total sensitivity indices. Hence the increase 
in computational cost is linear with the increase in the number of design variables. 

These indices can be used to understand the influence of design variables on the 
variability of any given objective over the chosen design space. This method proposed by 
Sobol (1993) is effective when the design variables are uncorrelated. In a multi-objective 
optimization study multiple optima are obtained which collectively give the Pareto 
optimal front (POF). The design variables at different regions along the POF share some 
similar features. Therefore, the design variables are correlated and this correlation has to 
be accounted for before estimating the influence of a design variable on an objective. The 
total sensitivity and main factor indices, which require uncorrelated data, cannot be used 
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for such a situation. Hence, for this purpose, a partial correlation coefficient (JMP, 2002) 
was calculated. 

Estimation of partial correlation coefficient (JMP, 2002) is a variance-based 
parametric approach that gives a measure of the linear relationship between the variances 
of a design variable, say xi and an objective J{X), after the influence of other variables 
have been filtered out. It gives a measure of expected change in_/(A) per unit change in xj 
which in other words gives the sensitivity of the function to the design variable. Linear 
approximations are obtained for*; and/(X) as a function of the components of Z and the 
residuals (say, ri and r 2 , respectively) measured as the difference between the data used 
for the approximations and the predicted value of the approximations are estimated. A 
partial correlation coefficient is the correlation between these two residuals and is given 
by 

fKvJiKvJi)] (3 . 27) 

<J cr 

r l r 2 

Optimization Algorithm 

The optimization process can be divided into two parts: 

1. RS A phase, 

2. Optimizer phase. 

In the first phase, RSA are generated with the available data set. In the second 
phase the optimizer uses the RSA during the gradient-based search for the optimum (the 
maximum or minimum of the objective). The initial design for the search is randomly 
selected from within the design space. The flowchart of the process is shown in Figure 
3.3. The optimization problem at hand can be formulated as min {/(X)} subject to lb < X < 
ub, where lb is the lower boundary vector and ub is the upper boundary vector of the 
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design variables vector X. If the goal is to maximize the objective function then/fx) can 
be written as -g(X), where g(X) is the objective function. Additional linear or nonlinear 
constraints can be incorporated if required. Solver (2002), an optimization tool available 
as part of Microsoft Excel package is an ideal tool for such studies. This tool uses the 
Generalized Reduced Gradient ( GRG2 ) nonlinear optimization code developed by 
Lasdon et al. (1978). 


r j 



Figure 3.3: Two phases of the optimization process, where Phase 1 deals with RSA and 
Phase 2 deals with optimization. 

In most optimization studies it is desirable to simultaneously optimize more than 
one response/objective. One method is to build, from the individual objectives, a 
composite objective known as the desirability function. The method allows for the 
designers’ priorities to be addressed during the optimization procedure. The first step is to 
define the desirability function, d for each objective. Desirability function is a 
normalizing of the objective such that its’ maximum and minimum over the design space 
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lie between zero and one. In the case where an objective should be maximized, say Ri, 
the desirability takes the form: 

d.= (3.28) 

\B-A) 

where B is the target value and A is the lowest acceptable value such that dj = 1 for any 
Rj> B and di = 0 for Rj < A. The weighting factor t is set according to the designer’s 
goal and the compromise he wants to achieve between different objectives. 

In the case where a response is to be minimized, say R 2 , the desirability takes on 
the form: 



(3.29) 


where C is the target value and E is the highest acceptable value such that d 2 = 1 for any 
R 2 <C and d 2 = 0 for R 2 > E. Choices for A, B, C, and E are made according to the 
designer’s priorities or simply as the maximum and minimum values of the objective 
over the design space. Values of s and t are set based on the priority of the objective. The 
sensitivity of the parameter, s, illustrated in Figure 3.4 can be instructive (Myers and 
Montgomery, 1995). Desirability functions with s « 1 imply that a product need not be 
close to the response target value, C, to be quite acceptable. However, s = 8, implies that 
the product is nearly unacceptable unless the response is close to C. A single composite 
desirability D is defined based on the geometric mean of the desirabilities of the 
individual objectives. 

D = {d r d 2 -d 3 ....d m f"> (3.30) 

This is then submitted to an optimization toolbox to be maximized. 
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Figure 3.4: Desirability function (d) for various weight factors, s. (Note: B < A) (Myers 
and Montgomery, 1995). 

Multi-Objective Genetic Algorithm 

One of the other optimization approaches is to use Multi-objective genetic 
algorithm (MOGA). It works on the principle of survival of the fittest. Genetic operators 
like reproduction, crossover and mutation are employed to find the optimal solution. In a 
multi-objective optimization scenario, when the objectives are conflicting in nature, many 
representative optimal solutions can be obtained. All these solutions comprise a set called 
the Pareto optimal set (Miettinen, 1999). The solutions in this set are considered non- 
dominated as they are not inferior to any other solution in the entire design space without 
having a bias towards some of the objectives. The function space of all the solutions in 
the Pareto optimal set is termed the Pareto optimal front (POF) (Miettinen, 1999). When 
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there are two objectives, the POF is a curve. When there are three objectives, the POF is 
represented by a surface. If there are more than three objectives, it is represented by 
hyper-surfaces. 

One of the MOGA that has been shown to be effective in finding the Pareto- 
optimal solutions is the elitist non-dominated sorting genetic algorithm (NSGA-II) (Deb 
et al., 2000). The algorithm can be described as: 

1 . Randomly initialize population (designs in the design space) of size N. 

2. Compute objectives for each design. 

3 . Rank the population using non-domination criteria. 

4. Compute crowding distance (this distance finds the relative closeness of the 
solution with other solutions in the objective space.) 

5. Employ genetic operators - selection, crossover & mutation - create new 
population. 

6. Evaluate objectives for the new population. 

7. Combine the two populations, rank them and compute the crowding distance. 

8. Select N best individuals. 

9. Go to step 3 and repeat till termination criteria is reached, which in the 
current study is chosen to be the number of generations. 

It is shown in a number of studies that using a combination of MOGA and local 
search (also known as hybrid GA), helps achieve faster convergence to the global Pareto 
optimal solution set (Deb and Goel, 2000, 2001; Goel, 2001; Goel and Deb, 2002; 
Ishibuchi and Yoshida, 2002). The a posteriori hybrid method (Deb and Goel, 2001) 
assumes the set of solutions generated by MOGA simulations as a starting point for the 
local search. Most of the local search methods are very efficient only for single objective 
optimization problems. Hence, a e-constraint strategy (Deb and Goel, 2001) helps reduce 
the dimensionality of the problem. Sequential quadratic programming, available in 
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Matlab (2002) can be used as the local search optimizer. This gives a set of solutions 
from which the dominated and duplicate solutions are removed to obtain the global 
Pareto optimal solution set and the POF. 

Hierarchical Clustering Method 

The POF can be divided into a number of clusters using a hierarchical clustering 
algorithm (Jain and Dubes, 1988) to assist the designer in selecting the optimal solution 
of choice. The clustering algorithm can be summarized as: 

1 . Start by assuming all the solutions as individual clusters. 

2. Find the mean of each cluster. 

3. Find the distance between clusters. 

4. Merge closest clusters. 

5. Go to step 2 till the number of clusters is equal to some prefixed value. 

6. Find the member of each cluster closest to the mean of the cluster. This is the 
representative element. 

Visualization Using Box Plot 

Box plot (JMP, 2002) is a visualization tool, which can help understand the 
variability in the design variables and objectives within the clusters of the POF. It can 
also assist in identifying the design variable that should be fixed when analyzing a given 
cluster. A schematic of a box plot with the y-axis representing the range of the scaled 
design variable is shown in Figure 3.5. 

Figure 3.5 shows that 25%, 50% and 75% of the data lie below the lower hinge, 
median and upper hinge, respectively. The difference between the upper and lower hinges 
is known as the “H-spread”. The inner fence is located 1 step beyond the hinges, which is 
equal to 1.5 times the H-spread. The upper adjacent value is identified as the largest value 
below the upper inner fence. The lower inner fence and lower adjacent value can be 
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similarly determined. These box plots will be used to visualize the variability along the 
POF. The spread between the upper and lower adjacent values gives the range of 
variation of a variable in a given set. Any variable lying beyond this spread is a potential 
outlier of the set. Small range gives a tight bound on that particular variable. Sometimes 
the box plot is known to collapse to a single point, which suggests that the variable 
should be fixed at that value. 



Figure 3.5: Schematic of a box plot 




CHAPTER 4 

DESIGN OPTIMIZATION - SINGLE ELEMENT ROCKET INJECTOR 

In this chapter the sensitivity and optimization studies of a single element liquid 
rocket injector design are presented. Initially a literature survey of the past injector design 
studies is presented. Following this the injector model is described. Then the results 
related to the sensitivity analyses and optimization studies are presented. 

Literature Review 

A critical goal for space propulsion design is to make the device safer, more 
affordable and more reliable. The design of combustion devices, namely, injectors, 
chambers and nozzles, will help in meeting these goals. The characteristics of the injector 
design are a key factor for both performance and thrust chamber environments. The thrust 
chamber performance is estimated by the rate and the extent to which mixing and 
resultant combustion occurs. The location of the mixing and resultant combustion 
determines the injector and thrust chamber thermal environments. These environments 
include temperatures on the combustor wall, the injector face and, for coaxial injectors, 
the oxidizer post tip. The difficulty encountered in designing injectors that perform well 
and have manageable environments is that the factors that promote performance often 
lead to increased heating of the solid surfaces of the injector and combustor thereby 
reducing the life expectancy or survivability of these components. 

Current injector design tools have been in use for 30 years or more and are largely 
empirical-based (Rupe, 1965; Dickerson et al., 1968; Pavli, 1979; Nurick, 1990; Pieper, 
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1991). Extensive sub- and full-scale hot-fire test programs often guided these injector 
design methodologies. The experimental databases, and thus the tools developed from 
them, are limited, in terms of design space, to specific element configurations that have 
been tested (Calhoon et al., 1973). In terms of scope, the design tools typically focus on 
performance, with the environment being a secondary consideration. The limited amount 
of environmental information available from these tools is usually one-dimensional and 
not functionally related to details of the injector design. It is very doubtful that 
application of these traditional design tools will result in robust future propulsion devices. 
Over time space propulsion programs with compressed schedules, lower budgets and 
more stringent requirements have resulted in the development of broader and more 
efficient injector design methodologies. 

The initial work by Tucker et al. (1998) and Vaidyanathan et al. (2000) focused on 
the optimization of a shear coaxial injector element with gaseous oxygen (GO 2 ) and 
gaseous hydrogen (GH 2 ) propellants. In this study, the design data was generated using 
an empirical design methodology developed by Calhoon et al. (1973). Calhoon et al. 
(1973) conducted a large number of cold-flow and hot-fire tests over a range of 
propellant mixture ratios, propellant velocity ratios and chamber pressure for shear 
coaxial, swirl coaxial, impinging, and premixed elements. The data were correlated 
directly with injector/chamber design parameters, which were recognized from both 
theoretical and empirical standpoints as the controlling variables. For the shear coaxial 
element, performance, as measured by energy release efficiency, ERE, was obtained 
using correlations taking into account combustor length, L CO mb (length from injector to 
throat) and the propellant velocity ratio, V/V a . The nominal chamber wall heat flux at a 
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point just downstream of the injector, Q„ 0 m, was calculated using a modified Bartz 
equation and was correlated with propellant mixture ratio, O/F, and propellant velocity 
ratio, V/V 0 to yield the actual chamber wall heat flux, Q. The objective was to maximize 
injector performance while minimizing chamber wall heat flux (lower heat fluxes reduce 
cooling requirements and increase chamber life/survivability) and chamber length 
(shorter chambers lower engine weight). 

In Tucker et al. (2000) the designs of an impinging injector element and a swirl co- 
axial injector element have been carried out. For the impinging injector element, the 
empirical design methodology of Calhoon et al. (1973) used the oxidizer pressure drop, 
AP 0 , fuel pressure drop, APf, combustor length, L com b, and the impingement half-angle, a 
as design variables. Objectives included ERE (a measure of element performance), wall 
heat flux, Q w , injector heat flux, Q inJ , relative combustor weight, W re i, and relative injector 
cost, C rd . The gaseous propellants were injected at a temperature of 540R. The empirical 
design methodology used to characterize the ERE and Q w used a quantity called the 
normalized injection momentum ratio to correlate the mixing at the different design 
points for the triplet element. 

The swirl coaxial element has been used sparingly in USA, but has been widely 
used in Russia because of its reported ability to perform well over a large throttle range 
(Gill and Nurick, 1976). This has sparked the interest in exploring its possibilities. The 
empirical design methodology of Calhoon et al. (1973) used the oxidizer pressure drop, 
AP 0 , fuel pressure drop, APf, , combustor length, L com b, and the full cone swirl angle, 0, as 
design variables. The objectives modeled were ERE (a measure of element performance), 
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wall heat flux, Q w , injector heat flux, Q in j, relative combustor weight, W re i, and relative 
injector cost, C re i- 

The advent and advancement of CFD in the last 20 years have produced a 
capability that has shown potential as an improved design tool in many areas of rocket 
propulsion. The application of CFD to injector design has lagged behind other areas such 
as turbomachinery because the physical models are more complicated for multiphase, 
turbulent reacting flows. New models that efficiently account for some of the complex 
processes (Cheng and Farmer, 2002) and thus increase the solution fidelity have recently 
become available. However, the three-dimensional geometry of multi-element injectors 
and the complex physical processes inherent in the flows that issue from them create 
major obstacles in validating the solution and generating significant amount of parametric 
solutions. The harsh high pressure and temperature environments typical of injector flows 
create significant difficulty in obtaining experimental data of satisfactory quality to 
validate and guide further development of computational models. Finally, solving the 
equations, for multiphase reacting flows, with high resolution typically requires lengthy 
computational times. However, continuing increases in computer speed and progress in 
parallel processing have begun to mitigate this turnaround problem. 

It has long been known that small changes in injector geometry can have significant 
impact on performance, as well as on environmental variables such as combustion 
chamber wall and injector face temperatures and heat fluxes (Gill and Nurick, 1976). 
Several geometric design variables can be accounted over appropriate ranges by 
calculating the performance and environmental variables of the injector using CFD. 
Global approximation method like RSM can provide a continuous representation of the 
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objectives over the design space thereby reducing the number of numerical computations 
required for sensitivity analyses. This global approximation can also guide the 
optimization study. Additionally, a global approximation method can also at times 
address issues related to CFD analysis, for e.g. grid resolution, thereby increasing the 
fidelity of the computations. The present work is first of its kind where a CFD-based 
design optimization methodology is proposed for the design of an injector (Vaidyanathan 
et al., 2003a; Vaidyanathan et al., 2004b). 

Scope 

The focus is on sensitivity and trade-off analyses for the design of a gaseous 
injector for liquid rocket propulsion. The data for the analyses is obtained from surrogate 
models (RSAs) of the objectives, and the multi-objective optima (Pareto optimal front, 
POF) generated by Goel et al. (2004) with the aid of multi-objective genetic algorithm 
(MOGA) and a local search method. The regions of the POF that represent different 
trade-offs among the objectives are obtained through a hierarchical clustering algorithm. 
The initial study concentrates on the generation of response surface approximations 
(RSAs) and preliminary optimization studies. Then it is followed by an elaborate 
sensitivity and trade-off analyses. The later part of the study is broadly divided into two 
parts. Firstly a sensitivity study is carried over the whole design space. The contribution 
of the design variables to the objective variability is calculated using a variance-based 
non-parametric approach (Sobol, 1993) and correlations between objectives are 
investigated. Secondly, sensitivity analyses are conducted on clusters along the POF. Box 
plots are used to highlight the variability of the design variables and the objectives within 
a given cluster. Additionally, the linear relationships between the variance of the design 
variables and the objectives are explored with the aid of partial correlation coefficients. 


67 


Injector model 

Liquid rocket propulsion injector elements can be categorized into two basic types 
based on propellant mixing. The first type is an impinging element (Figure 4.1 A) where 
mixing occurs by direct impingement of the propellant streams at an acute angle. The 
impingement enhances mixing by head-on interaction between the oxidizer and fuel (Gill 
and Nurick, 1976). The second type of injector consists of non-impinging elements where 
the propellant streams flow in parallel, typically in coaxial fashion (Figure 4. IB). Here, 
mixing is accomplished through a shear-mixing process (Calhoon et al., 1973). From a 
design standpoint, both element types have appealing characteristics. However both also 
have undesirable characteristics. For instance, if the impinging element has an F-O-F 
arrangement (Figure 4.1 A), the mixing occurs rapidly, which can yield high performance. 
However, since the combustion zone is close to the injector face, the potential for high 
levels of injector face heating must be considered. If the non-impinging element is 
assumed to be a shear coaxial element, mixing across the shear layer is relatively slow, 
requiring longer chambers to allow complete combustion. However, since the combustion 
zone is spread over a longer axial distance, the injector face is generally exposed to less 
severe thermal environments. 

Important design parameters for the impinging element (assuming fixed mass flow 
rates and constant propellant inlet conditions) include relative orifice size (or, relative 
stream momentum ratio), impingement angle and orifice spacing. Important parameters 
for the shear coaxial element (assuming fixed mass flow rates and constant propellant 
inlet conditions) include the area ratio of the two concentric tubes (or velocity ratio) and 
the shear area between the two propellant streams (i.e., the oxidizer post tip thickness) 
(Calhoon et al., 1973). 
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A B 


Figure 4.1: Schematic of the G02/GH2 impinging and coaxial injector elements. A) F-O- 
F impinging element. B) Coaxial element. 

One can combine the above-mentioned characteristics of injector types to develop 
hybrid concepts. For example, it has been noted (Gill and Nurick, 1976) that performance 
improvement in the shear coaxial element can be realized by directing the fuel toward the 
oxidizer stream rather than parallel to it and thinning the oxidizer post wall. The first 
modification causes the shear coaxial element to take on some of the aspects of the F-O-F 
impinging element. These notions lead to the hybrid element shown in Figure 4.2 that has 
been developed by The Boeing Company (U. S. Patent 6253539). 



Figure 4.2: Schematic of Hybrid Boeing Element (U. S. Patent 6253539) 


The four independent design variables chosen for the element used in this study are 


shown in Figure 4.3 A. They are the angle at which the H 2 is directed toward the oxidizer, 
a , the change in H 2 flow area from a baseline area, AHA, the change in 0 2 flow area from 
a baseline area, AO A, and the oxidizer post tip thickness, OPTT. The fuel and oxidizer 



















flow rates are held constant. The independent variable ranges considered in this exercise 


are shown in Table 4. 1 


0 2 Post Tip Thickness 
(OPTT=x-2x in) 


O 2 Flow Area 
(AOA=0~40%) 


H 2 Flow Area 
(AHA =0-25%) 


Figure 4.3: Design variables and objectives of the single element rocket injector. A) 
Design variables and their range. B) Objective functions. 
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Table 4.1: Range of design variables (a is an acute angle in degrees and x is the thickness 


of OPTT in inches) 

Variable 

Minimum 


Maximum 



Actual 

Scaled 

Actual 

Scaled 

a 

0° 

0 

a 0 

1 

AHA 

Baseline 

0 

Baseline+25% 

1 

AOA 

Baseline-40% 

0 

Baseline 

1 

OPTT 

x in 

0 

2x in 

1 


The dependent variables chosen for design evaluation are shown in Figure 4.3B. 
They are the maximum injector face temperature, TF max , the wall temperature at a 
distance three inches from the injector face, TW 4 , the maximum oxidizer post tip 
temperature, TT max and centerline axial location where the combustion is 99% complete, 
(hereafter referred to as the combustion length) Xc c (Figure 4.3B). The three temperatures 
(calculated as adiabatic wall temperatures in this study) were chosen as indicators of local 
environments. Lower temperatures would indicate a design that had longer 
life/survivability due to decreased thermal strain on the part. The combustion length, Xc C , 
was chosen as a measure of performance. Shorter combustion lengths indicate better 
performance. 

Numerical Procedure 

A pressure-based, finite difference, Navier-Stokes solver, FDNS500-CVS (Chen, 
1989; Wang and Chen, 1990; Chen and Farmer, 1991) is used in this study. The Navier- 
Stokes equations, the two-equation turbulence model, and kinetic equations are solved. 
Convection terms are discretized using either a second order upwind, third order upwind 
or a central difference scheme, with adaptively added second order and fourth order 
dissipation terms. For the viscous and source terms, second order central difference 
scheme is used. First order upwind scheme is used for scalar quantities, like turbulence 
kinetic energy and species mass fractions, to ensure positive values. Steady state is 
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assumed and an implicit Euler time marching scheme is used for computational 
efficiency. The chemical species continuity equations represent the H 2 - O 2 chemistry. It 
is represented with the aid of a 7-species and 9-reaction set (Chen, 1989; Wang and 
Chen, 1990; Chen and Farmer, 1991). The simulation domain and the boundary 
conditions used in all the CFD cases are shown in Figure 4.4. Because of the very large 
aspect ratio, both the injector and chamber have been shortened (at the cross hatched 
areas) for clarity. Both fuel and oxidizer flow in through the west boundary where the 
mass flow rate is fixed for both streams. The nozzle exit, at the east boundary, is modeled 
by an outlet boundary condition. The south boundary is modeled with the symmetry 
condition. All walls (both sides of the oxygen post, the outside of the fuel annulus, the 
outside chamber wall, and the faceplate) are modeled with the no slip adiabatic wall 
boundary condition. Each CFD analysis was done on a 200 CPU PC cluster with an 
AMD Athlon MP 1800 (1.8 GHz) chip and 1 GB RAM. All the cases were run 
concurrently and took about 5 days. 


adiabatic walk 


injector face 


oxygen po st 



'symmetry 

Figure 4.4: Simulation domain and boundary conditions 

Design Space 

In this study OA (with n c = 4, q = 2,t = 3 and X = 2) is used to generate 54 designs 


for RSA and testing of the approximation. To test the RSA, 14 of the 54 designs are 
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selected using cross-validation techniques (Papila, 2001). During the CFD simulations, 
two of the 40 designs (fitting set) were found to be unacceptable because they exhibited 
unsteady behaviors while the numerical algorithm used was a steady state model. Hence, 
the final information included 38 designs for fitting the RSA and 14 to test their 
predictive capabilities. All the design variables are scaled between zero and one based on 
their upper and lower bounds. All the objectives obtained from the CFD solutions of the 
52 valid designs are scaled between zero and one based on the maximum and minimum 
values and polynomial-based RS As generated. Once the RSAs are generated, the scaled 
objectives are normalized between zero and one based on the maximum and minimum of 
the generated RSAs, which will be referred to as the normalized objectives in the text. 
The scaled values are the non-normalized values of the objectives. The fitting and testing 
designs are shown in Tables 4.2 and 4.3, respectively. It should be noted that when the 
AOA is 1 or 0, the O 2 flow area is reduced by 0% or 40%, respectively, as compared to 
the baseline area. 

Results and Discussion 

CFD Analysis 

Comparison of two of the evaluated designs serves to illustrate the motivation for 
combining CFD analyses and an efficient global approximation model during the design 
process. The scaled independent design variables are shown for the two cases in Table 
4.4. In terms of the design space evaluated, these two designs are seen to be quite 
different. The normalized dependent variables are also shown in Table 4.4 with the 
temperatures shown in contour plots (Figures 4.5 and 4.6). 
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Table 4.2: Fitting designs (unacceptable designs are marked in bold). 


Case # 

a 

AHA 

AOA 

OPTT 

Case 1 : 

0 

0 

1 

0 

Case 2: 

0 

0 

1 

0.5 

Case 3: 

0 

0 

1 

1 

Case 4: 

0.5 

0.5 

0 

1 

Case 5: 

1 

1 

0.5 

0 

Case 6: 

1 

1 

0.5 

0.5 

Case 7: 

1 

1 

0.5 

1 

Case 8: 

0 

0.5 

0.5 

0 

Case 9: 

0.5 

1 

1 

0 

Case 10: 

0.5 

1 

1 

0.5 

Case 11: 

0.5 

1 

1 

1 

Case 12: 

1 

0 

0 

0 

Case 13: 

1 

0 

0 

0.5 

Case 14: 

1 

0 

0 

1 

Case 15: 

0 

1 

0 

0 

Case 16: 

0 

1 

0 

0.5 

Case 17: 

0 

1 

1 

1 

Case 18: 

0.5 

0 

0.5 

0 

Case 19: 

0.5 

0 

0.5 

1 

Case 20: 

1 

0.5 

1 

0 

Case 21: 

1 

0.5 

1 

0.5 

Case 22: 

1 

0.5 

1 

1 

Case 23: 

0 

1 

0.5 

0 

Case 24: 

0 

1 

0.5 

1 

Case 25: 

0.5 

0 

1 

0 

Case 26: 

0.5 

0 

1 

1 

Case 27: 

1 

0.5 

0 

0 

Case 28: 

1 

0.5 

0 

1 

Case 29: 

0 

0 

0 

0 

Case 30: 

0 

0 

0 

0.5 

Case 31: 

0 

0 

1 

1 

Case 32: 

0.5 

0.5 

0.5 

0.5 

Case 33: 

1 

1 

1 

0 

Case 34: 

1 

1 

1 

1 

Case 35: 

0 

0.5 

1 

0 

Case 36: 

0 

0.5 

1 

1 

Case 37: 

0.5 

1 

0 

0 

Case 38: 

0.5 

1 

0 

1 

Case 39: 

1 

0 

0.5 

0 

Case 40: 

1 

0 

0.5 

1 
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Table 4.3: Testing designs.. 


Case # 

a 


AOA 

OPTT 

Case 41: 

0.5 

0.5 

0 

0 

Case 42: 

0.5 

0.5 

0 

0.5 

Case 43: 

0 

0.5 

0.5 

0.5 

Case 44: 

0 

0.5 

0.5 

1 

Case 45: 

0.5 

0 

0.5 

0.5 

Case 46: 

0 

1 

0.5 

0.5 

Case 47: 

0.5 

0 

1 

0.5 

Case 48: 

1 

0.5 

0 

0.5 

Case 49: 

0.5 

0.5 

0.5 

0 

Case 50: 

0.5 

0.5 

0.5 

1 

Case 5 1 : 

1 

1 

1 

0.5 

Case 52: 

0 

0.5 

1 

0.5 

Case 53: 

0.5 

1 

0 

0.5 

Case 54: 

1 

0 

0.5 

0.5 


Table 4.4: Independent and dependent variable (objectives) for Cases X and Y from CFD 


computations (non-normalized values shown). 


Case 

a 

AHA 

AOA 

OPTT 

TF ma x 

tw 4 

TT m ax 

Xcc 

X 

1 

0 

0 

0 

0.998 

0.927 

0.128 

-0.004 

Y 

0 

0.5 

0.5 

1 

0.285 

0.395 

0.923 

0.567 



Figure 4.5: Temperature field for Cases X and Y. 


The chamber wall and injector face temperatures for Case Y (as seen in Figure 4.5) 
are low or moderately low, while for Case X, they are high. Figure 4.7 shows a large 
recirculation zone located between the injector and the chamber wall. This recirculating 
flow strips hot gases from the flame and causes them to flow back along the chamber 
wall and injector face. This phenomenon regulates the chamber wall temperature, TW4 
and the injector face temperature, TF max . Figure 4.6 shows that the other life/survivability 
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indicating variable, the maximum oxidizer post tip temperature, TT max , has essentially the 
opposite trend as compared to the other two temperatures. 



Figure 4.6: Near-injector temperature field for Cases X and Y. 



Large 

recirculation zone 


0 2 post tip 


wall 


Injector face 


Figure 4.7: Large recirculation zone in the combustion chamber. 
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The performance indicator, combustion length, X,* is seen (Table 4.4) to be at the 
minimum for Case X (shorter combustion lengths indicate better mixing elements) and at 
a moderate level for Case Y. Given these observations, it is clear that the dependent 
variables exhibit competing trends such that no design will produce the “best” values for 
all the dependent variables as desired in this study. 

These comparisons confirm the past observation that changes in the injector design 
details have major effects on injector performance and injector-generated environments. 
Injector designs which address space propulsion goals must be produced by a tool that 
accounts for performance and multiple, spatially resolved environmental variables. 
Efficient, validated CFD codes that model sufficient injector physics are necessary to 
meet these requirements. 

Generation of large amounts of complex information by these codes produces the 
need for a means to manage the data. The injector designer must be able to confidently 
and efficiently sort through this database to locate an acceptable compromise design. 
Global sensitivity and approximation tools can guide the designer effectively. 

Grid Sensitivity Investigation 

Initially the 54 cases identified by DOE were computed on an axi-symmetric 
geometry with 336x81 nodes. Only 33 out of the 40 fitting cases gave valid results. 
Results of the remaining seven cases contained unsteady features, which do not represent 
solutions of the steady state model employed. The RSAs generated with these data for 
TT max and X cc had R a 2 values of 0.961 and 0.976, respectively, suggesting a less than 
desirable fit. On checking the grid distribution in the combustion zone for the 33 cases 
used, it was determined that the grid resolution was insufficient. After a series of tests 
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involving addition of grid points in the axial direction in the combustion zone, a 430x81 
grid was found to be appropriate and used for the second run of the optimization study. 

To highlight the grid refinement, a comparison of the grid distributions is shown in 
Figure 4.8. The final grid was the product of tripling the axial node density in the 
combustion zone. The thick lines show the initial grid density, while the thin lines show 
final grid density. Note that, for clarity, only every sixth j-line is shown. New RSAs were 
generated from new solutions obtained on the fine grid. The new fits for TTma X and Xc C 
had R a 2 values of 0.989 and 0.995, respectively, representing a considerable improvement 
over the RS A performance based on the first, coarser grid. This time, only two out of the 
40 designs failed to produce valid results. It is to be noted that the R a 2 values for TF max 
and TW 4 are 0.999 for both the initial and final grid. This experience indicates that in 
addition to facilitating design optimization, the RSM can also help address the adequacy 
of the CFD solution accuracy. It offers insight into potential problems, based on the 
statistical regressions, from which the computations can be refined, thus improving the 
fidelity of the individual and collective databases. While this approach does not guarantee 
universally satisfactory outcomes, it does suggest clear directions to assess the critical 
area of data quality. 



Figure 4.8: Comparing the unrefined (336x81) {thick lines} and refined (430x81) {thin 
lines} grids. 
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Response Surface Generation 

In the following, data for each design variable from the 38 acceptable fitting cases 
were used to generate the RSAs. Both full and reduced quadratic polynomials are 
generated. Table 4.5 identifies o a and cr for the four objective functions and different 
polynomials generated. The full quadratic model is consistent in performance in terms of 
both <7 a and a. The reduced quadratic models either have a poorer value of <xor off er on jy 
marginal improvement over the full response surface model. Since there is no appreciable 
improvement by reducing the fits, there may be noise in the data or the quadratic models 
(both full and reduced) do not sufficiently represent the data. 


Table 4.5: Performance of full and reduced quadratic RSAs for the four objective 
functions (non-normalized values used). 




Full quadratic 

Reduced quadratic 

TF max 

Number of 
observations 

38 

38 ^ 


CT a 

0.00566 

0.00546 


a (14 points) 

0.00460 

0.00463 


Mean 

0.495 

0.495 

tw 4 

Number of 
observations 

38 

38 


<*a 

0.00803 

0.00795 


ct (14 points) 

0.00669 

0.00799 


Mean 

0.514 

0.514 

TT max 

Number of 
observations 

38 

38 


a a 

0.0413 

0.0401 


a (14 points) 

0.0396 

0.0382 


Mean 

0.560 

0.560 

Xce 

Number of 
observations 

38 

38 


CTa 

0.0205 

0.0197 


o (14 points) 

0.0178 

0.0186 


Mean 

0.497 

0.497 


Comparing the full quadratic model predictions for the fitting cases to the CFD 
results of the various objectives, the variations for TF max , TW 4 and Xc C were found to be 


79 


negligible (Figures 4.9A, 4.9B and 4.9D), suggesting no need for further improvement. 
However, a large number of points lie away from the best fit in the plot of TT max (Figure 
4.9C). 




Figure 4.9: Comparison between the best fit possible and as predicted by quadratic 

response surface. Optimum refers to the case when RSA and CFD values are 
identical. RS-CFD represents the value as for the current case (normalized 
Values Shown). A) TFmax- B) TW 4 . C) TTmax- D) Xcc- 

Based on this observation, a cubic model is generated for TT max . A full cubic fit in 
a 4 -design variable model requires a minimum of 35 designs points. To obtain a good fit, 
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the number of design points should be considerably larger than the required number. 
Since there are only 38 design points available from the fitting set, the testing set is also 
included in the fitting set. The PRESSrms is used to estimate the performance of the 
generated RSA along with the o a . Using the 52 design points, a full cubic was generated. 
The values of a a and PRESSrms were 0.0348 and 0.0598, respectively. The difference 
between the two measures of error is noticeable. Hence, a reduced cubic model is 
generated with the available design points. This improves the fit with values for cr a and 
PRESSrms being 0.0303 and 0.0388, respectively. Table 4.6 compares the performance of 
the full quadratic and the reduced cubic models for TT max . The reduced cubic model is 
seen to perform better statistically than the quadratic model. Hence, this cubic model is 
used for TT ma x in the optimization studies that follow. The RSAs for all four objectives 
are shown as Eqs. (4. l)-(4.4). These RSAs can then be used for sensitivity and 
optimization study. The RSAs are based on the scaled values of the design variables and 
objectives. 

Table 4.6: Performance of RSAs for the TT ma x- Reduced cubic RSA has 21 coefficients 


(non-normalized values used). 




Full quadratic 

Reduced cubic 

TT max 

Number of observations 

38 

52 



0.0413 

0.0303 


PRESSES 

0.0521 

0.0388 


Mean 

0.560 

0.591 


TFmax = 0.692 + 0.477(a) - 0.687(AHA) - 0.080(AOA) - 0.0650(OPTT) - 0. 167(a) 2 - 
0.0129(AHA)(a) + 0.0796(AHA) 2 - 0.0634(AOA)(a) - 0.0257(AOA)(AHA) + 


0.0877(AOA) 2 - 0.052 l(OPTT)(a) + 0.00156(OPTT)(AHA) + 0.00198(OPTT)(AOA) + 
0.0184(OPTT) 2 . (4.1) 
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TW 4 = 0.758 + 0.358(a) - 0.807(AHA) + 0.0925(AOA) - 0.0468(OPTT) - 0.172(a) 2 + 
0.0106(AHA)(a) + 0.0697(AHA) 2 - 0.146(AOA)(a) - 0.0416(AOA)(AHA) + 
0.102(AOA) 2 - 0.0694(OPTT)(a) - 0.00503(OPTT)(AHA) + 0.015 l(OPTT)(AOA) + 
0.0173(OPTT) 2 . (4.2) 

TT max = 0.370 - 0.205(a) + 0.0307(AHA) + 0.108(AOA) + 1.019(OPTT) - 0.135(a) 2 + 
0.0141(AHA)(a) + 0.0998(AHA) 2 + 0.208(AOA)(a) - 0.0301(AOA)(AHA) - 
0.226(AOA) 2 + 0.353(OPTT)(a) - 0.0497(OPTT)(AOA) - 0.423(OPTT) 2 + 
0.202(AHA)(a) 2 - 0.281 (AO A)(a) 2 - 0.342(AHA) 2 (a) - 0.245(AHA) 2 (AOA) + 

0.281 (AO A) 2 (AH A) - 0.184(OPTT) 2 (a) - 0.281(AHA)( a)(AOA) (4.3) 

Xcc = 0.153 - 0.322(a) + 0.396(AHA) + 0.424(AOA) + 0.0226(OPTT) + 0.175(a) 2 + 
0.0185(AHA)(a) - 0.0701(AHA) 2 - 0.251(AOA)(a) + 0.179(AOA)(AHA) + 
0.0150(AOA) 2 + 0.0134(OPTT)(a) + 0.0296(OPTT)(AHA) + 0.0752(OPTT)(AOA) + 
0.0192(OPTT) 2 (4.4) 

Optimization Process 

The generated RSA are used to conduct an optimization exercise and to study the 
relationship between the design variables and the objectives that are indicators of 
life/survivability and performance. Also, the ability to accommodate design features that 
promote extended life/survivability in a design while maintaining reasonable 
performance is explored. Three separate optimization studies, using the approach of 
desirability functions as explained in chapter 3, are presented below. First, single- 
objective minimizations are shown. Secondly, multi-objective optimizations are 
performed with equal weights. Finally, the multi-objective optimizations are conducted 
with variable weights. The obtained optimum solutions are compared with CFD 
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computations. Since the objectives are scaled to 0(1), the error measures have to be 
scaled accordingly to estimate the accuracy of the obtained solutions. In terms of the 
actual values the error for an objective is defined as 


error = 


_ ! yeFD }’RS 1 
y CFD 


(4.5) 


where y C FD is the solution obtained from the CFD analyses and )>rs is the prediction of the 
RSA. Using simple mathematics, not shown here, the error in the scaled variables can be 
written as 


error = 


y cfd y rs 

y CFD K 


(4.6) 


where the bar represents the scaled values, and K is defined as 


K = — ^ — (4.7) 

y max y min 

Here y m i„ and y max are the actual minimum and maximum values, respectively, 
based on the available set of fitting and testing data for that objective from the CFD 
analyses. In the dissertation, no bar is used for the notations used to represent the scaled 
design variables and objectives to avoid any confusion. Non-normalized values of the 
objectives refer to the scaled values and the normalized values refer to the re-scaled 
values of the objectives based on the maximum and minimum of the RSAs. This is to 
make some of the plots in the result section more meaningful. 

Single-objective optimization 

The purpose of this effort is twofold. First, the goal is to verify the performance of 
the optimization methodology in locating the minimum values for single objectives in the 
chosen design space. This verification is straightforward and necessary, but not sufficient 
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to conclude that the technique is useful for injector design. Table 4.7 shows the results of 
optimizing (here, minimizing) each objective separately. The numbers in parentheses in 
Table 4.7 indicate the weights applied during the optimization process. Here, a weight of 
one means that objective was included while a weight of zero means that objective was 
excluded. Based on the normalized values of the objectives the minimum value for each 


is necessarily 0. The optimizer found the correct value for each of the four cases. 


Table 4.7: Minimizing individual objectives. Value in parenthesis (1) indicates which 



Secondly, this process is helpful in understanding the injector operational 


influences via trend identification and variable groupings. The results from Table 4.7 are 


shown graphically in Figure 4. 10 to facilitate the discussion. It was shown earlier that the 
flow entrained in the large recirculation zone (Figure 4.7) regulated both the maximum 
temperature on the injector face, TF max and the chamber wall temperature, TW4. The 
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results from Opt-Cases 1 and 2 support this conclusion. Minimization of TF max very 
nearly minimizes TW 4 and vice versa. Not surprisingly, the designs for the two cases are 
also quite similar. Both designs are shear coaxial (a = 0) with the hydrogen flow area, 
AHA, and the oxidizer post tip thickness, OPTT, at their maximum values. Table 4.7 
shows that AOA is the only inconsistent variable in the two designs. In terms of the other 
objectives, when either TF max or TW4 is minimized, the maximum oxidizer post tip 
temperature, TT max , is high. The resulting moderate-to-long combustion lengths, Xc C , are 
consistent with the relatively slow mixing expected from a shear coaxial element. 
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Figure 4.10: Minimization of different objectives. (Case 1: TF max , Case 2: TW4, Case 3: 
TT max , Case 4: Xc C . Normalized values shown). A) Objectives. B) Design 
variables. 

Reference to Opt-Cases 3 and 4 in Table 4.7 and Figure 4.10 also suggests a 
correlation between TT max and Xc C , although not as tight as the one for TF^x and TW4. 
When TT max is minimized, Xc C is low and vice versa. Further investigation needs to be 
done to completely understand the physics that underlies this correlation. Here, both 
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designs are impinging-like with the hydrogen flow angle, a, at or near its maximum 
value and AHA and OPTT at their minimum values in the chosen design space. Again, 
Table 4.7 shows that AOA is the only inconsistent variable between the two cases. Both 
minimizations result in very high values for TF max and TW 4 . Not surprisingly, the 
impinging elements represented by the two designs yield significantly shorter combustion 
lengths than the shear coaxial designs. 

The trends for a, AHA and OPTT are consistent between the two pairs of 
dependent variables but AOA varies among the four cases. With the other design 
variables set at the optimum design levels, Figure 4.1 1 A shows the variation of TF max and 
TW 4 with AOA. A similar plot for TT max and Xc C is shown in Figure 4.1 IB. The 
maximum injector face temperature, TF max , is least sensitive to AOA. It is also the only 
objective to exhibit a minimum value in the interior of the design space. This finding 
suggests that expanding the design space could result in more robust designs for the 
multi-objective optimization. The trends for TW 4 and Xc C are similar to each other, with 
AOA zero for both cases. The trend for TT max relative to AOA is opposite. 

Table 4.7 and Figure 4. 10 show the resulting design is different for each of the four 
single-objective optimizations. These results indicate that designing robust injectors 
demands consideration of all the independent variables together. Figure 4.10 clearly 
shows the competing design trends. Thus, meeting a set of design requirements will 
require a compromise design. 



86 



Figure 4.1 1: Variation of objectives with respect to AOA for fixed a, AHA and OPtt 

(normalized values shown and D is the desirability function. All objectives are 
equally weighed in the composite function). A) TF max and TW 4 for a= 0 , 
AHA=1 and OPTT=l. B) TT max and X« for <x=l, AHA=0 and OPTT=0. 

Multi-objective optimization 

To concurrently evaluate component life/survivability and performance 
considerations, a multi-objective optimization study is carried out. Starting with 
performance, X cc , as most of the injector design tools do, the objectives influencing 
thermal environments are added one at a time to study the effect on the resulting 
optimum designs. These optimum designs are presented in Table 4.8 and Figure 4.12 
where Opt-Case 4, with Xc C minimized is repeated as the starting point. The design for 
this case is an impinging element (a = 0.917) with minimum flow areas and the thinnest 
oxidizer post tip. The consequences of this design are a minimum Xc C , a low TT max and 
very high values of TF max and TW4. 

When TT m ax is minimized along with X cc in Opt-Case 5, the optimum design (foes 
not change significantly from that in Opt-Case 4. The hydrogen flow angle increases 
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form 0.917 to 1.0. The values of AHA, AOA and OPTT remain unchanged from Opt- 
Case 4. The combustion length, X cc , is unchanged as compared to Opt-Case 4, whereas 
TT max is marginally improved from 0.181 to 0.152. The other thermal objectives remain 
at very high levels. Recall from Table 4.7 that AOA was 1.0 for the single-objective 


optimization of TTmax. 

Table 4.8: Study of effect of life/survivability on performance. Value in parenthesis 



a 

AHA 

AOA 





Xcc 
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0.917 

0 



mu 

mtm 

0.182 
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CFD 
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0.981 

eeqh 


-0.004 
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0.12 

5 

1 

0 

0 

0 

1.0 
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mm 

mi 

CFD 





0.998 

0.927 


-0.004 






0.03 

0.01 


0.14 

6 

1 

1 

0 

0 

0.376 

(1) 

0.224 
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0) 

mm 
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0.214 

0.252 

0.264 









0.38 

i 

1 

1 
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mm 
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0.369 

0.214 

0.252 

0.264 

Error 
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0.10 

0.12 

3.93 

0.38 


More insight regarding the optimum design can be gained by visualizing Figure 
4.1 IB. Although the individual TT max and X cc drive the design towards the opposite ends 
of the range of AOA, the composite desirability function, D, has marginal variation with 
respect to AOA. The optimizer picks the value of AOA that gives the maximum value of 
the desirability function, D. The optimum design at this point is affected mostly by the 







































88 


other design variables. This shows the benefits of using a global optimization procedure, 
which helps identify the trend of the composite function over the design space even when 
the individual objectives have opposing trends. 
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Figure 4.12: Composite minimization of objectives with different weightings. (Case 4: 
(0,0,0, 1), Case 5: (0,0, 1,1), Case 6: (1,0, 1,1), Case 7: (1,1, 1,1). The values in 
parenthesis indicate weights for (TF max , TW 4 , TT max , X cc ). Normalized values 
shown). A) Objectives. B) Design variables. 

Opt-Case 6 minimizes TF max simultaneously with TT max and X cc . In terms of the 
design variables, AHA has now shifted from 0 to 1, while the other three remain at the 
same values as in Opt-Case 5. This shift in AHA is consistent with the single-objective 
minimization of TF max shown as Opt-Case 1 in Table 4.7. This single change in the 
design produces dramatic changes in the objectives. As desired, the value of TF max has 
decreased significantly from 1.0 to 0.376. The concurrent large drop in the value of TW 4 
is consistent with the earlier conclusion linking TFmax and TW 4 . The oxygen post tip 
temperature, TT max , remains essentially unchanged. The increase in X cc from 0 to 0.279 is 
a negative impact of the design change. Figure 4.13 sheds more light on the situation. 
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Figure 4. 13 A plots the individual objectives as a function of AHA. Figure 4.13B is 
similar, with the individual objectives plotted as a function of AOA. These figures clearly 
illustrate that both TF max and TW 4 are considerably more sensitive to AHA than AOA. 
Accordingly, the optimizer uses the hydrogen flow area to most efficiently regulate 
TF max . Figure 4.13A also shows the positive slope of Xc C relative to AHA, which is 
responsible for the fairly large performance drop seen in this design. 




Figure 4. 13: Variation of objectives with respect to AHA or AOA with other design 

variables fixed, (normalized values shown and D is the desirability function. 
All objectives are equally weighed in the composite function). A) TF max , TW4, 
TT max and Xc C , with respect to AHA for a=l, AOA=0 and OPTT=0. B) TF max , 
TW4, TT max and Xc C , with respect to AOA for a=l, AHA=1 and OPTT=0 

Opt-Case 7 adds the final objective, TW4. Table 4.8 shows that addition of TW4 has 
no effect on the design and consequently no effect on any of the objectives. The chamber 
wall temperature, TW 4 , has previously been shown to be linked with TF raax , so this result 
is not surprising. 

The injector design resulting from inclusion of all the objectives in the optimization 
process (Opt-Case 7) is different from Opt-Case 4 where only the performance indicator, 
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Xcc was optimized. In the design space evaluated for this effort, AHA has a very strong 
effect on the system, especially for TF max and TW4. Reference to Figure 4.13 shows that 
at the zero value for both AHA and AO A, X cc is still decreasing. Performance could 
probably be further improved (i.e., Xc C shortened) by enlarging the design space to 
include smaller values of AHA and AOA. In the process of designing a real injector, the 
large drop in performance affected by the consideration of TF max and TW4 would 
probably require a compromise design with certain variables weighted more heavily than 
others. 

To further probe the interplay between injector performance and component 
life/survivability, a composite optimization study is conducted with varying weights on 
the individual objective functions embodied by the composite desirability function. This 
is a functionality that the new design optimization methodology must possess. The 
designer must be able to weight, or favor, one or more dependent variables to have 
maximum flexibility in addressing design requirements. The use of geometric mean and 
desirability functions is an approach which addresses such a requirement. Additionally 
MOGA have also been shown to work well for such problems. The results of the MOGA 
will be presented in the later sections of this chapter. The results of the current study are 
shown in Table 4.9 and Figure 4.14. Reference to Table 4.9 (the numbers in parenthesis 
indicate the weight in the composite desirability function for that variable) shows that the 
designs suggested by Opt-Case 8 and Opt-Case 10 are similar. In Opt-Case 8 the 
life/survivability influencing objectives TF max , TW 4 , and TT max are weighted by a 2:1 
ratio over the performance indicator (X^). Opt-Case 10 is the opposite, with performance 
having a 2:1 weighting over the life/survivability objectives. Reference to Tables 4.8 and 
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4.9 shows that the design variables and objectives for these two cases are almost identical 
to those in Opt-Case 7, where all four objectives are weighted equally. Figure 4. 13B 
shows the composite desirability function D, for Opt-Case 7, to be reasonably flat as a 
function of AOA. Thus, small weightings have little or no effect on the injector design. 
Table 4.9: Study of influence of life/survivability and performance objectives on each 


other. Value in parenthesis indicates the weighting on the objective functions 
(normalized values shown). 
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ii 

0.612 

0 

0 




0.282 

(0.1) 

EH 

CFD 





0.911 

0.890 

0.361 

0.0070 

!■ 





0.10 

0.09 

3.10 

0.16 


The other two cases shown in Table 4.9, Opt-Case 9 and Opt-Case 11, have relative 
weightings of 50:1 and 1:50, respectively for life/survivability (TF raax , TW 4 , and TT max ) 
and performance indicators (X cc ). For Opt-Case 9, the design tends toward the cases 
where TF ma x and TW 4 were individually minimized. The exception is that OPTT is equal 
to zero, which is required to minimize TT max . For Opt-Case 1 1, the resulting design tends 
toward Opt-Case 4 where Xc C is minimized. Even with this strong weighting, the effect of 
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including TF raa x and TW4 is felt with the value of a equal to 0.612. When X ec is 


minimized by itself, a is equal to 0.917. 
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Figure 4.14: Composite minimization of objectives with different weightings. (Case 8: 
(1,1, 1,0.5), Case 9: (5, 5, 5, 0.1), Case 10: (0.5, 0.5, 0.5,1), Case 11: 

(0.1,0. 1,0. 1,5). The values in parenthesis indicate weights for (TF max , TW4, 
TT max , Xc C . normalized values shown). A) Objectives. B) Design variables. 

Figure 4.15 shows the joint desirability function, D, for Opt-Case 9 (with a 
weighting of 50: 1 in favor of life/survivability over performance), to be fairly flat. 
Accordingly, small weightings have little effect on the design. Significant weightings 
must be applied to push the design very far toward either life/survivability or 
performance. Opt-Case 1 1 has good performance and a low injector tip temperature. 
However, the chamber wall and injector face temperatures are very high. In an actual 
design, if this level of performance was required, active cooling could be required for the 
chamber wall and injector face. 

Additional CFD solutions were executed to confirm the optimum designs obtained 
from RSA for all 1 1 Opt-Cases. These results are shown in Tables 4.7, 4.8 and 4.9. The 
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objective error, as defined by Eq. (4.5), is also shown for each case. The error is on the 
order of four percent or less in the case of TT mBX . The error for the other three objectives 
is less than one percent. The larger discrepancy in TT max is likely due to the steady state 
assumption made in the CFD analyses for this effort. The injector problem is actually 
unsteady. 


Response Functions 



Figure 4.15: Variation of TF max , TW 4 , TT max and Xc C , with respect to AOA for a=0, 
AHA=1 and OPTT=0 (normalized values shown and D is the desirability 
function. Case 9: Temperatures weighted over performance in the ratio of 
50:1). 

Preliminary Observations 

Based on the preliminary optimization studies specific observations can be made. 
Single-objective optimization. The optimizer reliably locates single objective 
minimum values. Minimization of the individual objectives yields a different injector 
design. This establishes the fact that designing an efficient long life/survivability injector 
requires concurrent evaluation of all the relevant design variables. 

For the current study, 
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1 . The four objectives, based on the design obtained during their individual 
minimization, largely fall into two groups: (TF max , TW 4 ) and (TT max , Xc C ). 
However, there are differences between them. In each group, the responses of 
the two design objectives to AOA are different. The observation also 
indicates that the three life/survivability-related objectives require 
compromises between design variables. 

2. Minimizing TF max and TW 4 leads to a design with a equal to zero (shear 
coaxial element), maximum fuel flow area and thickest post tip. This design 
also yields moderate to poor performance due to the slow mixing across the 
shear layer. 

3. Minimizing TT max and X cc results in an impinging-like design with a equal to 
one. It also has the minimum fuel flow area and the thinnest post tip 
thickness. This design performs well, but has very high wall and injector face 
temperatures. 

Multi-objective optimization. The injector design, when multiple objectives are 
included, is different from any of the single objective minimizations. For example, 
although the individual objectives, TT max and X cc drive the design towards the opposite 
ends of the range of AOA, the composite desirability function accounting for all their 
effects exhibits only marginal variation with respect to it. The optimum design is 
affected mostly by other design variables. 

One can clearly see the benefit of using a global optimization procedure, which 
helps identify and interpret the trend of the composite function over the design space, 
especially when the individual objectives have opposing patterns. 

For the current study the following specific observations can be made. 

• Equal weights 

1 . The design with all 4 objectives included (Opt-Case 7) is still an impinging 
element with low to moderate temperatures, but marginal performance. 

2. Figure 4. 13 indicates that lowering either propellant flow area beyond the 
current limit would probably decrease Xc C and thus increase performance. 
Lowering the oxidizer area would probably have less adverse effect on TF max 
and TW 4 . 
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• Unequal weights 

1 . With modest weights on either the performance (Xc c ) or life/survivability 
(TFnmx, TW 4 and TT max ) the design does not change appreciably. 

2. Opt-Case 1 1 (with a large emphasis on performance) gives very good 
performance, modest TT max , but very high values of TT max and TW 4 . If 
enlarging the design space does not help, active cooling may be required. 

Sensitivity and Trade-Off Analysis 

The analysis is divided into two parts. Firstly, the global analysis is carried out 
where the entire design space is explored to estimate the sensitivity of the objectives to 
the design variables. Following this the POF is explored to understand the trade-offs 
between the objectives and also the different design trends. 

Global analysis 

The global analyses addresses, both the sensitivity of objective variability to design 
variables and the interaction between objectives over the complete design space. The 
global sensitivity indices are computed using the variance-based non-parametric 
approach proposed by Sobol (1993) and described in chapter 3. The code written using 
Matlab (2002) was validated using a well-known benchmark problem (Mckay, 1997). 
Table 4. 10 summarizes the results of this study and lists the essential and non-essential 
design variables with respect to individual objective variability. A design variable is 
considered essential if it is responsible for at least 5% of the objective variability. 

Figure 4.16 shows the percentage of main factor (Si) contribution of different 
design variables to individual objective variabilities. The variability of TF max is largely 
influenced by AHA and moderately by a (Figure 4.16A). The effect of the other design 
variables is marginal suggesting that they are non-essential and in principle could be 
fixed. The variability of TW 4 is considerably influenced only by AHA (Figure 4.16B). 
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TT max is influenced considerably by OPTT and marginally by a (Figure 4. 16C). For X cc , 
AHA, AOA and a have considerable influence (Figure 4.16D). The total sensitivity 
indices ( S ‘ otal ) were computed and compared with the main factors (Si) and it was found 
that the contributions of the cross-interactions among the design variables to the 
objectives variability were negligible. 

Table 4.10: List of essential ( V) and non-essential (x) variables for each objective and the 


mean errors between the modified RSA and original RSA. An essential 
variable accounts for at least 5% of the objective’s variability. 


Objective 


Design 

a 

Variables 

AHA 

AOA 

OPTT 

Mean 
Error (%) 

Life/Survivability 

TF max 

V 

V 

X 

X 

5.2 

indicators 

tw 4 

X 

V 

X 

X 

11.5 


TTmax 

V 

X 

X 

V 

6.7 

Performance 

indicator 

Xc C 

V 

V 

V 

X 

6.1 


TFmax TW4 


OPTT, 



AHA, 0.955 


A B 

Figure 4.16: Main factor (Si,) influence on objective variability. A) TF max . B) TW 4 . C) 
TT m ax- D) X cc . 
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TTmax 


Xcc 



Figure 4.16: Continued. 

This global sensitivity analysis shows that different design variables can have 
varied effect on individual objectives. Such a study can help the designer fix some of the 
design variables while analyzing the design. For example, in the current injector design 
study, for objective TF max , the design variables AOA and OPTT can be fixed at their 
mean value (0.5) as this do not result in significant differences in the prediction. The 
mean error (difference between predictions) throughout the design space between 
modified RSAs (obtained by fixing the non-essential design variables at their mean 
values), and the RSAs listed in Eqs.(4.1)-(4.4), (referred to as original RSA in the rest of 
the text) are given in Table 4.10. The error for TW 4 is about 12% suggesting that the cut- 
off for the non-essential variables may have to be lowered to capture additional features 
of the original model. The mean error in the modified RSA of the remaining objectives is 
about 6%, suggesting that they capture the original RSA reasonably well. This 
information can be used for dimensionality reduction and therefore, to ease the search for 
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optimum designs. Similar concept has been addressed in the past using a different 
approach. For example, Knill et al. (1999) have used linear aerodynamic to identify the 
important terms in the polynomial-based RSA which were then used for creating the 
surrogate model from Euler analyses. This reduced considerably the number of points 
needed for the Euler design of experiments. For the current work since the number of 
design variables is small the remaining studies are carried out without fixing the non- 
essential variables. From an engineer’s viewpoint and interest, this study provides an 
insight into the physics of a design problem by highlighting the design features that 
govern the individual objectives. 

Linear surrogate models for the four objectives are constructed as a function of the 
four design variables. The corresponding coefficients are shown in Table 4. 1 1 . The 
magnitude of the coefficients agrees well with the essential and non-essential nature of 
the design variables for each objective. Additionally, the R 2 values for the linear RSA are 
compared with those of the original RSAs (referred as R 2 nonlinear in Table 4.1 1). 
Comparing R 2 unear with R 2 nonlinear shows that most of the variability is accounted for by 
the linear model and the additional terms in the nonlinear model give marginal 
improvement. 

Table 4.1 1 : Coefficient associated with the different terms in the linear RSA and 


comparison ofi? 2 values of linear RSA with that of nonlinear RSA. 



a 

AHA 

AOA 

OPTT 

Intercept 

R linear 

p2 

nonlinear 

TFmax 

0.237 

-0.634 

-0.033 

-0.066 

0.733 

0.990 

0.999 

tw 4 

0.0719 

-0.775 

0.114 

-0.052 

0.826 

0.986 

0.999 

TT max 

-0.222 

0.108 

-0.063 

0.662 

0.367 

0.918 

0.994 

Xcc 

-0.234 

-0.402 

0.433 

0.108 

0.138 

0.958 

0.997 
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A correlation analysis was then carried out by Goel et al. (2004) to observe the 
interactions between objectives. A total of 14641 (1 1 4 ) design points were generated over 
the complete design space by varying one variable at a time by a constant value. The 
objectives were calculated using the original RSAs. The correlation matrix Cdes and 
corresponding p- values are computed using Matlab. 

The correlation matrix, Cdes shows that there is a strong correlation between 
objectives TF max and TW4, as the corresponding coefficient is very close to 1. P-values 
and 95% confidence intervals for the correlation coefficients also establish the statistical 
significance of the results. Low P-values (« 0.05) confirms the significance of the 
correlation results. This finding is in agreement with the observations made in the 
preliminary optimization study. The combustion chamber wall temperature, TW4, is 
excluded and the optimization problem is formulated with the remaining three objectives. 

' TF max X cc TW 4 TT^ 

TF max 1.00 -0.773 0.947 -0.272 
Cdes = x„ -0.773 1.00 -0.583 0.263 
TW 4 0.947 -0.583 1.00 -0.193 
TT -0.272 0.263 -0.193 1.00 

max 

Pareto front analyses 

The trade-offs between objectives and sensitivity analyses are carried out at the 

POF. The Pareto fronts of (TF max , Xc C ) and (TF max , TTmax) are of interest as these 
objectives conflict one another. Goel et al. (2004) have generated the POFs with the aid 
of NSGA-II. Figure 4.17A shows the relation between TFmax and Xc C . It can be seen that 
the POF in this case is linear over a large region. A small increase in the value of TF max 
(« 10%) reduces the combustion length, X cc , by nearly 50%. Figure 4.17B shows the 
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relation between TF max and TT max . It is obvious that the POF is non-convex. It is also 
seen that a small drop in the value of the face temperature (« 10%) can reduce the tip 
temperature TT max by nearly 60%. Hence at a small cost of TF max both X cc and TT max 
improves considerably. 



A B 


Figure 4.17: Two-objective Pareto optimal front. A) TF max vs Xc C . B) TF max vs TT max . 

Following the trade-off studies, the three-objective (TFmax, Xc C and TT max ) 
optimization problem is solved using the NSGA-II algorithm (Goel et al., 2004). The 
solutions obtained are further refined using the e-constraint local search strategy coded 
in Matlab. In this strategy, first, one of the objectives TF max , is treated as the objective 
function and rest of the objectives X« c and TT max are treated as equality constraints. The 
constraint value is set equal to the corresponding objective value as found by NSGA-II 
simulation. The corresponding design variable vector is used as initial guess. This 
procedure is repeated for all individuals in the population. This gives a set of Pareto 
optimal solutions, referred to as Set A in this study. Next, Xc C is chosen as the objective 
and the other objective functions TT max and TF max are treated as constraints. For this 
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problem, the obtained Pareto optimal set is referred to as Set B. Similarly, the Pareto 
optimal set, Set C, is obtained with TT max as objective and TF max and X cc as constraints. 
Now all the Pareto optimal solutions obtained so far; i.e.. Sets A, B, C and the original 
NSGA-II set are combined. To find the true Pareto optimal front, a non-domination check 
is carried out on this set of 400 solutions. This yields 254 non-dominated solutions. After 
removing the duplicates, there are 249 solutions, which are on the POF. These solutions 
are shown in Figure 4.18. This solution set is the global Pareto optimal solution set. The 
optimal solutions listed in Tables 4.7-4.9 are also plotted in Figure 4.18 



Figure 4.18: Pareto optimal solution set (+) and the multi-objective optima listed in 
Tables 4.7-4.9 (o). 

The optimal solutions shown in Tables 4. 7-4.9 lie along the POF obtained using 
MOGA and the local search algorithm. Although identical solutions are not identified by 
the different multi-objective optimization approaches, the trend of the Pareto front is 
adequately captured. This shows the effectiveness of the different approaches. 
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Goel et al. (2004) has then used the hierarchical clustering algorithm (Jain and 
Dubes, 1988) to divide the obtained POF into 9 clusters for sensitivity and trade-off 
analyses. Values of the design variables and objectives for these designs are shown in 
Table 4.12. Graphically, the solutions are shown on the POF in Figure 4.19. It is clear 
from Figure 4.19 that solutions are selected uniformly over the design space. In Figure 
4.20, the box plots of the design variables for clusters 1, 3, 6 and 9 are shown. The box 
plots for the objectives are shown in Figure 4.21. These plots highlight the variability of 
the design variables and objectives in each cluster. 


Table 4.12: Objective function and design variables of nine (9) representative designs 


from the Pareto optimal solution set. 

Cluster 

a 

AHA 

AOA 

OPTT 

TFmax 

Xc C 

TTmax 

1 

0.000 

1.000 

0.842 

0.712 

0.023 

MEM 

0.880 

2 

0.000 

1.000 

0.356 

0.587 

0.028 

Mm 

0.0890 

3 

0.000 

1.000 

0.442 

0.014 

0.054 

wttm 

0.452 

4 

0.094 

1.000 

0.000 

0.015 

0.126 

0.453 

0.466 

5 

0.668 

1.000 

0.732 

0.000 

0.259 

0.681 

0.229 

6 

0.600 

0.670 

0.000 

0.000 

0.489 

0.264 

0.226 

7 

0.295 

0.108 

0.000 

0.354 

0.719 

0.129 

0.641 

8 

0.314 

0.066 

0.000 

0.055 

0.776 

0.097 

0.357 

9 

1.000 

0.014 

0.680 

0.000 

0.935 

0.138 

-0.043 


For cluster 1 it is seen that the value of a is fixed at 0 (shear coaxial injector) 
(Figure 4.20A) and AHA is fixed at 1 (Figure 4.20B). This suggests that in this cluster, 
the designs are sensitive to a and AHA both of which reach their extreme values. The 
remaining two design variables AOA (Figure 4.20C) and OPTT (Figure 4.20D), vary 
over a range. It is observed that TF ma x is minimized (Figure 4.21 A) where as Xc C and 
TT m ax lie near their maximum and have little variability. Hence the designs in cluster 1 
tend to minimize TF max and represent shear coaxial injector designs. The design variables 
AOA and OPTT do not influence TF max but affect the remaining objectives, Xc C and 
TT max . Partial correlation coefficients are estimated to obtain the relationship between the 
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variances of these design variables and objectives (Table 4.13). It is noticed that as AO A 
increases, Xcc increases (Rcoit = 1-00) and TT max decreases (Rcon = -0.638), thereby 
requiring a compromise in the design. As OPTT decreases both Xc C and TT max decrease 
(Rcoit = 1 -00 for both). 



Figure 4.19: Pareto optimal solution set and nine (9) representative solutions from the 
same set. The circles identify the representative of a specific cluster. 




Figure 4.20: Box plots for the design variables in clusters 1, 3, 6 and 9. A) a. B) AHA. C) 
AOA. D) OPTT. 
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Figure 4.20: Continued. 

Similar observations can be made for clusters 3, 6 and 9. Figures 4.21A-4.21C 
show that with increase in TF max , from cluster 1 to cluster 9, Xcc and TT max decrease 
(based on the median of the box plots). This highlights the trade-off between the 
objectives. Cluster 9 provides information about the opposing trend as to what was 
observed in cluster 1. As objectives Xc C and TT max are minimized (Figures 4.21B and 
4.2 1C) an impinging injector design is obtained (a ~ 1, Figure 4.20 A) with TF max 
exhibiting high values (Figure 4.21 A). The AHA is near minimum (Figure 4.20B) 
contrary to the design in cluster 1 where high AHA minimized TF max . The AOA has 
considerable variation which suggests that the variability of objectives, X cc and TT max are 
not largely affected by this design variable. Figure 4.20D shows that Xc C and TT max are 
minimized for the minimum value of OPTT. 

Table 4.13 gives the partial correlation coefficients for the set of design variables 
and the objectives in each cluster which shows considerable variation. The partial 
correlation coefficients for the combinations left out are effectively zero. These measures 
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can be used to tune the design variables in a chosen cluster so as to improve on the 


objectives as per design requirements. 





Figure 4.21: Box plots for the objectives in clusters 1, 3, 6 and 9. A) TF max . 


TT 


max- 


B)Xc C . C) 


In this study an elaborate optimization study with different compromised design 


goals has been presented along with both a global and POF sensitivity analyses. The 
conclusion drawn from these studies will be presented in chapter 5. 
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Table 4.13: Partial correlation coefficients (Rcorr) of design variables vs. objectives for 


different clusters along the POF 

Cluster 

Design variable 

TF m ax 

x cc 

TT m ax 

1 

AOA 

- 

1.000 

-0.638 

1 

OPTT 

- 

1.000 

0.991 

3 

AOA 

- 

1.000 

- 

6 

a 

0.982 

-0.735 

-0.983 

6 

AHA 

-0.999 

0.994 

- 0.729 

9 

a 

0.877 

-0.203 

- 0.769 

9 

AHA 

-0.992 

0.983 

0.816 

9 

AOA 

-0.977 

0.997 

-0.940 


CHAPTER 5 

SUMMARY, CONCLUSION AND FUTURE WORK 

In this chapter the work done is this dissertation is summarized and conclusions 
drawn. Firstly, the summary and conclusion of the work related to Navier-Stokes (NS) 
CFD code verification is presented. Following this the design optimization study is 
summarized and conclusions drawn. Finally possible future research topics are briefly 
mentioned. 

Navier-Stokes Code Verification 

In this study, least square extrapolation (LSE) method has been mainly 
implemented on Navier-Stokes computations solved in the finite volume (FV) 
formulation. Two coarse grid solutions are extrapolated onto a fine grid using constant 
weighting parameters estimated by minimizing the residuals of the NS equations on the 
fine grid in the least square sense. The lid-driven cavity flow with two different boundary 
conditions (constant and variable lid velocities) are used to study the robustness and 
efficiency of LSE in situations of singularities, non-linearity and coupled equations. The 
initial test for LSE is done using a 2-D turning point problem is tested. Following this, 
only pressure is extrapolated for the lid-driven cavity flows. Finally LSE is implemented 
on the complete NS equations to extrapolate the velocity components and pressure 
simultaneously. 

The key observations pertaining to the LSE approach can be summarized as 
follows: 
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• It has been noticed with the aid of a 2-D turning point problem and a laminar lid- 
driven cavity flow that this method can yield solution improvement for linear as 
well as nonlinear PDEs. Additionally it has been shown with this method that the 
extrapolating weights based on minimization of the residual of the governing 
equation performs better than the weights used by RE based on assumed order of 
convergence. 

• It has been shown with the aid of two different Reynolds number (Re = 5.33 and 
1000) that LSE works well over a wide range provided the base grids captures the 
flow features adequately. 

• Minimizing the residuals does not always lead to minimizing the extrapolated 
solution error. In particular, it has been shown that flow with singularities Would 
introduce large variation in the eigenvalues of the set of residual equation that is 
used for least square approach. This would influence the accuracy of the solution 
even though least square might effectively reduce the L 2 norm of the residual. 

• In practice, neglecting the regions of singularity can improve the performance. This 
suggests the requirement of further study to address issues of singularities in 
different flows. 

• The LSE is seen to work well for the pressure-velocity coupled system. The Picard- 
like iteration adapted for the nonlinear momentum equation converges well f or the 
shown case. 

• It is also seen that the individual extrapolation of pressure or the coupled pressure- 
velocity extrapolation does not affect the accuracy of the extrapolated pressure field 
drastically. 

Overall issues related to singularities in the lid-driven cavity flow were noticed. 
This resulted in an elaborate study to understand the problems related to it and its 
influence on LSE. In this study only a constant weighting parameter is used. The 
preliminary implementation of spatially dependent weighting function exhibits 
inconsistent performance, which requires further investigation. 

Design Optimization Study 

In this study, the integration of CFD and optimization tools to explore the design of 
a single element rocket injector is presented. The design variables indicating the thermal 
environment and performance of the element are identified and the design points selected 
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using a DOE technique. CFD solutions are obtained for each design and appropriate 
orders of polynomial-based RSAs are generated for the individual objectives. Based on 
the performance of the RSAs a grid refinement study is carried out and a new, more 
efficient grid is generated. Utilizing the results obtained from CFD computations on the 
fine grid, a preliminary multi-objective optimization study is carried out using a 
composite function based on geometric mean approach and different design trends, 
correlations and trade-offs observed. 

Following the preliminary study, RSA, global sensitivity indices and genetic 
algorithms have been used to conduct an elaborate sensitivity and trade-off analyses of 
the injector design. Global analyses have been conducted to observe the influence of 
design variables on the objective variability and estimate the correlation between 
objectives. A Pareto optimal front (POF), generated by Goel et al. (2004) using a RSA- 
multi-objective genetic algorithm (MOGA) coupled approach, is used to conduct trade- 
off studies between objectives. A three-objective POF has been divided into 9 clusters 
and box plots have been used to observe the variability of the designs in each cluster. 
Additionally, partial correlation coefficients have been calculated to derive a relationship 
between the objectives and design variables in each cluster. Such analyses provide the 
designer with enough information to filter out designs based on his specific needs. Some 
observations based on the tools used, their performance, and the different design aspects 
of the injector are highlighted below. 

• The degree of spread of a dependent variable over a parametric space is an 
important criterion for validation experiment design. Evaluation of large parametric 
spaces provides insight into potential validation measurements. 

• Correlations of the variables highlight the trends and provide insight into physical 
processes that regulate the objectives. 
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• The study confirms that the details of the injector design govern the performance. 
Directing the fuel flow towards the element axis and thinning the oxidizer post lead 
to better performance. The same details are shown to have a large impact on the 
injector and combustion chamber environmental variables as well. 

• Using CFD to evaluate and include more design variables and objectives in the 
conceptual design phase has potential to create more robust designs. This increase 
in the dimensionality of the design problem also creates the requirement for an 
optimization tool to guide the designer through often opposing trends to an 
acceptable design compromise. 

• Based on the preliminary optimization study it was found that the proposed CFD- 
RSM approach worked well. It was fairly efficient in terms of the number of CFD 
solutions required and provided sensible results in view of the physics of the 
problem. 

• For both single- and multi-objective optimization efforts, in many cases, the 
resultant designs are at an extreme of the design space. This indicates the design 
space could be enlarged subject to practical design considerations. Again, a merit 
of the RSM is that the solutions obtained can be used repeatedly when revising the 
design scope. Furthermore, different scenarios, with either single- or multi- 
objective optimization, can make use of the information supplied by the response 
surface repeatedly. 

• The global sensitivity indices (main factor and total sensitivity indices) help 
identify which design variable is essential for objective variability. Additionally, 
point out the effect of cross interactions among the design variables on the 
objectives variability. 

• The variability of TF max depends on a and AHA, TW4 on AHA, TT max on a and 
OPTT and Xc C on a, AHA and AOA. The rest of the design variables for each 
objective can, in principle, be fixed, as their contribution to the objective variability 
is marginal. It is noticed that the mean errors between the modified RSAs and the 
original RSAs are about 6-12%. 

• It is observed that TF max and TW 4 are strongly correlated and hence TW 4 is 
excluded from the multi-objective optimization study. 

• The conflicting nature of the objectives (TF max , TT max ) and (TF max , Xc C ) is observed 
by examining the corresponding POFs. For a small cost in TF max (~10%) , both Xc C 
and TT max (~50%) can be decreased considerably. 

• The POF obtained for the 3-objective study (TFmax, TT max and Xc C ) gives an idea of 
the various compromises among the different objectives. A hierarchical clustering 
approach helps divide the POF into regions of similar objective goals. 
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• Box plots and partial correlation coefficients can assist the design selection along 
the POF. Box plots identify the variability of design variables and objectives in 
each cluster. Partial correlation coefficient identifies linear trends between the 
design variables and objectives variabilities. 

• To minimize TF ma x it is essential to have a shear coaxial injector (a = 0) and 
maximum change in baseline H 2 flow area. To minimize Xc C and TTmax, an 
impinging injector (a > 0) is required with minimum change in baseline H 2 flow 
area. 

T his study has offered in depth information about the different design aspects 
pertaining to the design variables and objectives. The importance of different design 
variables on individual objectives has been identified and the interaction between the 
objectives noticed. The POF of the 3-objective study gives different choices for the 
designer to look into but the key observation is that for a marginal cost in the injector 
face temperature, TF max , the performance, X cc and O 2 post tip temperature, TT max can be 
considerably reduced. 

The MDO tools, namely DOE and RSM, have been shown to work well along with 
CFD. Vaidyanathan et al. (2003b) have shown that these tools can be use to address 
model fidelity issues in the context of a turbulent cavitating flows. In this work they have 
used these tools to identify the optimum values of the empirical parameters present in the 
turbulence and the cavitation models. This highlights the broad scope of the tools 
presented in this dissertation. 

Future work 

Future research directions are briefly presented in this section. First issues related 
to LSE are addressed and possible strategies suggested. Next, ideas for effective use of 
design tools used in this dissertation are presented. 
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Navier-Stokes Code Verification 

In this dissertation, the weight function in LSE has not been allowed to vary 
spatially. For the problems tested and the techniques implemented, spatially dependent 
weight function was found to deteriorate the extrapolation and hence only a constant 
weight was used for all the study presented. One of the key issues during the 
extrapolation is the interpolation of coarse grid solutions onto the fine grid. This 
interpolation introduces high spurious wave number terms which combined with the 
spatially dependent weight function deteriorates LSE. For flows like lid-driven cavity 
flows, such high wave number components are amplified during residual estimation. The 
LSE approach is implemented on these polluted residuals and hence does not guarantee 
the minimization of the solution error. More details regarding this issue can be found in 
the work of Garbey and Shyy (2004). They have proposed using a low mode 
approximation of the residual for LSE implementation. This smoothing of the residual is 
shown to be effective in improving the performance of LSE. 

One of the other issues of LSE implementation is that the coarse grids should 
capture the flow feature adequately for efficient extrapolation. At the same time these 
coarse grids should not be so fine as to increase the computational cost. As the grid 
resolution is increased the solution reaches an asymptotic trend; i.e., there is no 
noticeable change in the flow feature with further refinement of the grid. Unless the 
coarse grids are properly chosen, the performance of LSE will not be optimal. Judicious 
selection of the base grids is an issue that deserves to be investigated. 

Design Optimization 

In this dissertation the RSAs were first generated and then used for the sensitivity 
and trade-off analyses. Hence although the global sensitivity analyses provided incite into 
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the influence of design variables on the individual objective variability, it was not used 
for dimensionality reduction. Additionally the injector design problem has small number 
of design variables and hence there is not a serious requirement for such a reduction. For 
future design studies the following direction is suggested to tackle the curse of 
dimensionality. Whenever possible the preliminary design study can be conducted using 
a reduced order model, including semi-empirical and experimental data. This reduced 
model can be coupled with DOE, RSM and global sensitivity tools to identify the 
influence of design variables on objective variability. This information can then be used 
to reduce the dimensionality of the design space. Higher fidelity models can then be used 
to analyze the designs in this reduced design space. Additionally, the reduced design 
space will help in reducing the number of search directions during the optimization study. 
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